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Research  Overview 

This  report  serves  as  the  final  summary  of  our  research  under  the  Netflow  program  as 
part  of  the  Composite  CAD  Program.  In  this  report  we  will  discuss  the  proposed 
objectives  of  the  Netflow  program  and  the  accomplishments  over  the  course  of  the 
program. 

2.  Program  Summary 

The  NetFlow  program  is  aimed  at  developing  Computer-Aided  Design  (CAD)  tools  for 
the  analysis  of  microchemical  fluidic  networks.  The  specific  aim  of  the  program  is  to 
develop  software  tools  for  the  analysis  of  the  diverse  transport  phenomena  that  are  present 
in  these  systems  -  these  include  pressure  driven,  electrokinetic  and  diffusive  transport. 
Experimental  validation  for  the  software  tools  developed  is  an  inherent  part  of  this 
program. 

The  milestones  and  deliverables  for  the  program  are  outlined  in  the  Gantt  Chart  in 
Table  1.  Variances  from  the  originally  proposed  milestones  are  also  indicated  in  the  chart. 

The  overall  program  summary  is  as  follows: 

•  All  milestones  originally  set  forth  in  the  program  have  been  met.  Some  of  the  tasks 
were  completed  well  ahead  of  schedule,  whereas  some  others  were  delayed  for 
various  reasons,  but  we  were  able  to  accommodate  the  delays  by  faster  progress  in 
other  areas.  Details  of  each  task  and  results  are  attached  in  the  appendix. 

•  Two  additional  tasks  were  added  in  the  third  year,  beyond  the  scope  of  the  original 
program  -  the  first  is  tools  for  the  analysis  of  non-dilute  sample  transport,  and  the 
second  tools  for  the  analysis  of  pneumatic  plug  driven  flows  which  is  an  additional 
transport  mechanism  in  microchemical  networks.  These  tasks  were  also  completed 
within  the  time  frame  specified. 

•  The  NetFlow  software  has  been  released  commercially.  Three  versions  of  the 
software  have  been  released  in  the  past  two  years  -  each  version  is  an  enhancement 
of  the  previous  one.  Tools  are  available  for  the  analysis  of  containment  (NetFlow- 
C),  transport  (NetFlow-T)  and  for  the  design  of  microchemical  networks  (NetFlow- 
G).  The  software  is  available  on  Sun,  HP,  SGI  and  Windows  NT  Platforms.  The 
User’s  Guide  and  Reference  Manuals  give  much  more  details  about  the  software 
and  its  use.  These  are  attached  in  the  appendix. 

•  The  software  allows  the  analysis  of  mixed  transport  mechanisms,  including 
pressure-driven,  electrokinetic,  pneumatic  and  diffusive  transport.  More  complex 
transport  mechanisms  including  the  non-dilute  sample  transport  is  also  incorporated 
in  the  software. 

•  The  software  tools  are  in  continued  development  with  additional  releases  planned 
for  the  next  several  years.  Current  users  include  both  industrial  and  academic 
groups. 
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•  NetFlow-G  is  currently  available  commercially  and  allows  generator  driven  creation 
of  shapes  for  device  layout. 

•  A  caged  fluorescence  experimental  setup  was  developed  at  PE  Applied  Biosystems 
and  data  was  gathered  for  validating  the  developed  tools.  Stanford  also  collected 
data  at  Sandia  using  the  caged  fluorescence  setup  there.  Stanford  has  undertaken  the 
development  of  their  own  set  up  in  the  past  year  and  the  first  results  are  now 
flowing  from  their  setup. 

•  NetFlow  has  found  application  in  helping  the  design  of  current  and  next-generation 
devices  and  device  components  in  PE- Applied  Biosystems  and  Stanford,  as  well  as 
in  other  research  groups. 

•  Data  analysis  tools  were  developed  in  the  program  to  enable  the  validation  of 
simulation  and  experimental  results.  These  tools  are  available  as  part  of  the 
software. 

•  Several  research  publications  and  presentations  have  emerged  over  the  course  of 
this  program.  These  include  papers  at  Hilton  Head  ‘98,  Transducers  99  and  the  |i 
conferences,  among  others.  A  complete  list,  including  copies  of  the  papers  is 
attached  in  the  appendix. 

3.  Program  milestones 
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NetFlow 


Introduction 

The  NetFlow  module  aids  the  user  in  the  simulation  and  design  of  interconnects  for  chemical  transport  of  species  in 
microchannels.  The  physical  phenomena  associated  with  the  transport  include: 

■  Pressure  flow 

■  Diffusion 

■  Electrophoresis 

■  Electroosmosis 

NetFlow  can  be  applied  to  model  these  phenomena,  either  individually  or  in  combination.  The  equations  used  to 
model  the  flow  are  the  Navier-Stokes  equations  for  the  flow,  coupled  with  the  mass  equation  for  the  species  motion 
and  the  Poisson  equation  for  the  electric  field.  Electrophoretic  effects  are  modeled  through  the  application  of  an  elec¬ 
trophoretic  mobility  that  represents  the  force  on  the  charged  species  due  to  the  electric  field.  Electroosmotic  effects 
can  be  modeled  by  the  direct  solution  of  the  Poisson-Boltzmann  equation  for  the  electrical  double  layer,  or  alterna¬ 
tively  through  the  specification  of  an  electroosmotic  mobility  coefficient  for  the  walls. 

The  NetFlow  module  is  very  useful  for  designing  the  “wires”  for  chemical  transport,  for  designing  the  containment  of 
a  sample  in  processing  nodes,  and  for  making  detailed  quantitative  measurements  of  microchemical  flows  to  calibrate 
and  validate  simulations. 
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NetFlow  Theory 

Introduction  to  NetFlow 

The  NetFlow  module  is  built  on  the  base  MemCFD  module  in  MEMCAD  and  is  specific  to  the  analysis  of  microflu¬ 
idic  transport  channels  in  chemical  and  biological  analysis  systems.  On-chip  analysis  systems  consist  of  active  and 
passive  components  as  shown  in  the  figure.  Examples  of  the  former  are  reaction  chambers,  switching  joints  etc., 
whereas  the  latter  include  transport  and  separation  channels. 


For  a  particular  system  of  interest  the  design  analysis  generally  involves  the  representation  of  the  device  as  a  system 
comprising  several  individual  components.  The  physics  of  each  component  is  then  extracted  from  detailed  computa¬ 
tion  for  that  component.  Each  individual  component  can  thus  be  represented  in  the  system,  yielding  a  reduced  order 
model  for  the  device.  This  approach  is  analogous  to  the  circuit  simulation  methodology  adopted  in  both  mechanical 
and  electrical  domains. 

Fluidic  systems,  however,  differ  in  one  significant  aspect  from  mechanical  and  electrical  systems.  In  both  mechanical 
and  electrical  systems  the  interconnections  between  the  components  can  be  considered  “ideal”  - — that  is,  the  output  of 
one  component  can  be  directly  applied  as  the  input  to  the  next.  In  the  fluidic  domain,  this  is  not  generally  true.  The 
interconnects  here  are  channels  that  transport  the  chemicals  from  one  processing  site  to  another,  often  under  the  effect 
of  diffusion  and  other  related  phenomena.  As  a  result,  interconnects  often  have  a  profound  effect  on  the  behavior  of 
the  chemical  (the  “signal”  in  the  system).  Analyzing  these  effects  requires  modeling  of  the  channels. 

In  summary,  the  design  problem  is  as  follows:  For  a  given  application,  the  various  physical  properties  that  affect  the 
behavior  of  the  device  must  be  extracted  in  an  appropriate  parameter  space.  The  properties  can  then  be  used  in  that 
parameter  space  to  conduct  further  analyses.  NetFlow  is  the  first  step  in  enabling  this  process. 
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Governing  Equations 

The  basic  equations  describing  the  fluid  motion  are  the  Navier-Stokes  equations  with  appropriate  electromigratory 
flux  terms  to  represent  the  effect  of  the  applied  electric  field  on  the  carrier  and/or  the  charged  species.  The  effect  of 
the  applied  field  can  be  divided  into  two  fundamental  components,  electrophoresis  and  electroosmosis. 

Eiectrophoresis 

The  basis  for  electrophoresis  is  the  differential  migration  of  the  charged  species  ions  relative  to  the  carrier  molecules 
under  the  application  of  the  external  field.  The  differential  migration  is  primarily  an  effect  of  the  difference  in  the  net 
charge  between  the  solvent  and  solute  ions,  although  frictional  effects  may  also  have  some  relevance.  The  migration 
velocity  of  the  charged  species  can  be  expressed  in  terms  of  the  applied  field  strength  as 

Vep  =  v-epE  m 

where  |Tep  is  the  electrophoretic  mobility  of  the  ion  in  the  carrier  species.  It  is  important  to  note  that  in  most  cases 
the  carrier  does  not  move  under  electrophoresis. 

Electroosmosis 

Electroosmosis,  in  contrast,  is  a  macroscopic  phenomenon  involving  the  pumping  of  a  fluid  through  a  channel  under 
the  application  of  the  field.  In  most  cases,  walls  in  microchannels  are  characterized  by  the  presence  of  surface 
charges.  The  charge  may  either  be  due  to  the  property  of  the  wall  or  by  adsorption  of  the  charged  species  from  the 
buffer.  In  the  presence  of  an  electrolyte,  the  surface  charge  density  induces  the  formation  of  a  double  layer  in  the  fluid 
by  attracting  oppositely  charged  ions  from  the  electrolyte  to  the  immediate  vicinity  of  the  wall.  The  application  of  the 
electric  field  exerts  a  force  on  the  fluid  that  is  initially  felt  only  within  the  double  layer.  As  a  result,  the  fluid  in  the 
near  vicinity  of  the  wall  starts  to  move.  Due  to  the  viscous  forces,  the  fluid  in  the  center  of  the  channel  is  also  acceler¬ 
ated  until  the  net  velocity  gradient  in  the  radial  direction  is  zero  and  the  whole  fluid  in  the  channel  moves  at  a  con¬ 
stant  velocity. 

The  determination  of  the  electroosmotic  flow  field  requires  the  solution  of  the  Navier-Stokes  equation.  Incorporating 
the  electroosmotic  effect  as  a  force  the  equations  become 

v  •  v  =0  (Continuity) 
eo 

D(p  V) 

- —  =  -VP  +  V-aVV  +oE  [2] 

Dt  eo  e 


Where  Veo  is  the  induced  electroosmotic  velocity,  D/Dt  is  the  material  derivative  of  the  momentum,  |T  is  the  fluid 
viscosity  and  P  is  the  pressure.  The  last  term  in  the  momentum  equation  represents  the  electroosmotic  force  on  the 
fluid.  Here  pe  is  the  charge  density  and  E  is  the  electric  field  intensity.  The  electric  field  can  be  determined  by  the 
solution  of  the  potential  equation 

E  =  -VO  V20  =  —  (Pe  /  £)  [3] 
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The  electric  potential  O  can  be  further  decomposed  into  two  components,  (j),  the  external  applied  potential  and  the 
zeta  potential  at  the  walls. 


0  =  4)  +  ^  i4] 


Under  the  assumption  that  the  zeta  potential  effects  are  confined  to  a  very  small  region  near  the  wall,  the  charge  dis¬ 
tribution  can  be  assumed  to  be  governed  by  the  zeta  potential  alone,  independent  of  the  external  field.  This  allows  the 
decomposition  of  Eq.  [3]  into  separate  equations  for  the  applied  field  and  the  zeta  potential 


V >  =  0  (Applied  potential)  [5a] 

2 

V  C,  —  — (p^/s)  (Zeta  potential)  [5b] 


The  equation  for  the  zeta  potential  can  be  simplified  by  the  Debye-Hiickel  treatment  for  the  charge  density.  This  is 
based  on  the  concept  of  a  diffuse  double  layer,  proposed  initially  by  Gouy  and  Chapman,  which  assumes  that  the  dou¬ 
ble  layer  extends  for  some  finite  distance  into  the  fluid.  Through  the  Debye-Hiickel  approximation  the  charge  density 
on  the  walls  can  be  determined  through  the  Boltzmann  equation  and  results  in  the  Poisson-Boltzmann  equation  for 
the  zeta  potential 


[6] 


where  I1q  is  the  ionic  concentration,  Z  is  the  valence  of  the  charged  buffer,  Q  is  the  electron  charge,  k  is  the  Boltz¬ 
mann  constant  and  T  is  the  temperature.  Eq  is  the  permittivity  of  vacuum. 

Under  further  simplification  under  the  assumption  that  ZQ  «  kT  equation  [6]  reduces  to  the  form 


2  2 

VC  =  k  C 


where 


K  = 


(2  «0z 


2  2 
e 


/  E ZqUT) 


1 

2 


is 


known  as  the  Debye-Hiickel  parameter  and  1/K  is  termed  as  the  characteristic  thickness  of  the  double  la 


ayer. 


The  coupled  solution  of  the  Eqs.  [2],  [5a]  and  [7]  yields  the  flow  field  under  the  conditions  specified. 
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The  motion  of  a  charged  species  in  the  electric  field  can  be  determined  by  coupling  the  species  mass  equation  [8]  with 
the  above  equations. 


dt 


+  (Veo  +  VVC*  =  °^Ci 


[8] 


The  convective  transport  of  the  species  is  through  the  combined  effect  of  the  electroosmotic  motion  of  the  carrier 
fluid  and  the  electrophoretic  transport  of  the  species  under  the  effect  of  the  applied  electric  field. 


The  equations  presented  in  this  section  are  derived  under  the  following  assumptions: 

■  Dilute  Carrier:  The  carrier  fluid  is  assumed  to  be  electroneutral  everywhere,  except  within  the  double 
layer.  The  charge  distribution  is  therefore  confined  to  a  small  region  within  the  double  layer. 

■  Dilute  Solution:  The  carrier  fluid  is  assumed  to  be  the  predominant  species  in  the  mixture.  That  is,  the 
mass  fraction  of  the  charged  species  is  assumed  negligible  in  comparison  to  the  fraction  of  carrier.  This 
is  generally  valid  in  most  electrokinetic  problems  where  the  species  under  observation  are  generally  in 
the  millimolar  range  (corresponding  to  mass  fractions  of  1CT3  and  lower).  This  implication  of  the  dilute 
assumption  is  that  the  solute  species  does  not  affect  the  material  properties  of  the  mixture.  This  allows 
the  species  transport  equation  to  be  decoupled  from  the  momentum  equation. 

■  Individual  Species  do  not  affect  each  other:  This  assumption  is  used  in  describing  the  flux  terms  in 
the  above  equations,  in  that  the  flux  of  one  species  does  not  depend  on  the  fluxes  of  the  other  species. 
This  is  only  true  in  dilute  solutions. 

■  No  Chemical  Reactions:  The  charged  species  are  assumed  to  be  fully  ionized  in  the  mixture. 

As  mentioned  above,  the  dilute  solution  allows  the  density  of  the  mixture  to  be  assumed  constant  and  equal  to  the 
density  of  the  carrier.  This  reduces  the  problem  to  the  incompressible  form.  The  decoupling  of  the  species  equation 
from  the  momentum  equation  also  implies  that  the  species  conservation  condition  is  not  relevant  and  can  be  dis¬ 
carded. 
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Appendix  1:  NetFIow  Users  Guide 
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Reference 

NetFlow  window  map 
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NetFlow  -  Model  window 


The  Model  window  displays  the  current  mbif  files  and  scratch  workspace;  the  window  is  the  same  as  viewed  with 
other  MEMCAD  4  functions.  It  is  not  updated  automatically  from  other  windows. 
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□  MBIF  File 

The  base  model  file  used  for  the  solver  calculations.  By  convention,  the  solver  always  uses  the  base  model 
mbif,  even  though  mbif  files  created  by  NetFlow  and  other  solvers  contain  all  the  information  needed  to  com¬ 
pute  a  solution.  The  path  is  not  set  automatically  from  other  screens. 

□  Scratch  workspace 

The  path  used  for  storing  temporary  files  created  during  the  solution.  Any  temporary  files  may  be  deleted 
when  the  solver  has  completed  its  task. 

□  Run  Name 

The  path  used  to  store  the  single  result  mbif  file  created  by  NetFlow.  The  directory  name  should  provide  some 
way  of  identifying  the  type  of  run  created. 

□  OK 

Close  the  window  and  return  control  to  the  MEMCAD  4  main  menu. 
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NetFlow  -  Tool  Setup  window 

The  Tool  window  allows  user  to  set  the  parameters  for  running  the  fluidic  solver.  Most  of  the  information  in  the  tool 
window  provides  the  “how”  and  “what”  instructions  to  NetFlow. 


□  Controls/Method 


The  Navier-Stokes  equations  used  in  the  FlumeCAD  solvers  are  non-linear  and  require  an  iterative  approach 
to  solving  problems.  The  solvers  adjust  a  variety  of  variables  as  part  of  the  solution  method,  allowing  them  to 
solve  both  straightforward  problems  with  single  solutions,  as  well  as  more  complex  problems  with  multiple 
local  minimum  or  maximum  values. 

■  Coupled:  The  coupled  solver  is  the  default  choice  for  FlumeCAD.  This  iterative  solver  is 
recommended  for  most  problems,  and  can  handle  steady  state  and  transient  flows. 

■  Segregated:  The  segregated  solver  solves  each  degree  of  freedom  in  the  problem  separately  and 
consecutively.  Thus,  the  three  momentum  equations  for  flow  in  the  X,  Y,  and  Z  directions  and  the 
pressure  and  energy  equations  are  each  solved  separately.  This  approach  uses  less  system  memory  at  the 
expense  of  a  longer  solution  time.  It  is  recommended  for  problems  containing  more  than  2500  elements. 

□  Controls/Advanced  Controls 

Several  different  tolerances  may  be  set  depending  on  the  non-linearity  of  the  flow  problem. 

■  Solution  Tolerance:  Convergence  value  setting.  For  large  problems,  the  default  value  of  0.01  ensures  a 
solution  in  a  reasonable  time.  For  smaller  problems  this  value  may  be  decreased  even  further. 

■  Residual  Tolerance:  Similar  to  the  Solution  Tolerance;  controls  the  Flux  of  a  particular  quantity. 

■  Upwinding:  A  numerical  solution  stabilization  technique  that  prevents  solver  instability.  Occasionally, 
the  solver  will  be  unable  to  complete  its  calculations  due  to  mesh  anomalies.  If  this  occurs,  change  this 
value  to  1  (second-order  stabilization)  or  2  (first-order  stabilization). 

■  Relaxation:  Specifies  relaxation  factors  for  the  degrees  of  freedom  in  the  system.  The  default  window 
setting  of  0.0  maintains  a  0.05  factor  for  all  degrees  of  freedom.  The  number  should  be  changed  only  to 
modify  the  degree  of  coupling  between  the  system  equations.  Large  factors  of  0.2-0.4  can  be  used  for 
strong  coupling,  while  lower  numbers  are  best  for  weak  coupling.  Optimal  values  are  problem- 
dependent,  and  can  be  derived  only  by  trial  and  error. 

■  Memory  Use:  The  default  setting  is  - 1 ,  which  allows  MEMCAD  4  to  allocate  memory  as  required.  For 
large  3D  problems,  or  to  avoid  unnecessary  swapping,  the  user  should  set  this  value  if  the  job  is  much 
bigger  than  the  RAM  memory.  A  suggested  range  is  2/3  to  4/5  of  actual  RAM. 

□  Transients/Analysis 

■  Steady:  Performs  a  steady  state,  or  time-independent,  analysis.  The  solver  will  try  and  converge  to  an 
answer  within  the  number  of  steps  given  in  the  Edit  Steady  window  Maxlncrement  setting. 

■  Transient:  This  setting  sets  up  the  transient  simulation.  Five  time  values  are  set  in  the  Edit  Transient 
window: 

Start:  Time  to  start  simulation.  This  value  is  usually  0. 

Delta_T:  Time  interval  for  simulation.  This  is  the  initial  timestep  when  TimeStep  is  set  to  Variable, 
and  is  the  actual  timestep  when  TimeStep  is  set  to  Fixed.  In  either  case,  it  should  be  set  to  a  very 
small  number  (roughly  le-6). 

Stop:  Time  to  stop  simulation. 

DTMax:  When  TimeStep  is  set  to  Variable,  sets  the  maximum  range  for  the  time  step.  If  set  to  0,  it 
is  ignored.  It  should  be  set  to  0,  except  for  complex  problems,  where  it  should  be  set  to  1/100  of 
the  total  simulation  time  (0.01  *  (Stop  -  Start)). 

Output:  Specify  time  increments  at  which  the  mbif  files  are  written.  For  example,  a  setting  of  0.2 
writes  an  mbif  file  every  0.2  seconds  for  the  duration  of  the  simulation. 
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□  Transients/TimeStep 

This  controls  the  time  step  of  the  overall  transient  analysis. 

■  Variable:  This  MEMCAD  4  time  stepping  control  default  value  is  the  optimal  setting.  The  solver 
determines  the  time  step  value  at  any  given  point  in  the  solution.  With  a  variable  setting,  the  solver  can 
use  small  time  steps  at  the  beginning  of  the  problem,  and  larger  steps  near  the  end  when  the  solution  is 
nearly  complete. 

■  Fixed:  This  time  step  is  used  for  special  problems.  An  example  is  when  the  solver  would  increase  the 
timestep  in  the  variable  mode  with  steps  that  are  larger  than  the  time  in  which  the  effects  that  are  studied 
occur.  In  highly  viscous  flows  the  actual  startup  effects  of  the  flow  might  otherwise  be  skipped  with  a 
variable  setting.  However,  note  that  compute  times  may  increase  dramatically.  As  an  illustration,  if  the 
Start  time  and  Stop  time  in  the  Transient  field  is  set  to  0  and  1  second,  and  the  Delta_T  to  le-6,  then  in 
theory  it  would  take  the  solver  le6  steps  in  order  to  calculate  the  answer.  This  setting  works  well  with 
the  Coupled  Method  setting. 

□  Transients/Tolerance 

Tolerance  of  a  physical  time  step  in  the  transient  analysis.  Depending  on  the  mass  conservation  plot  results, 
the  Tolerance  value  may  be  increased  or  decreased  from  the  default  setting  of  0.001 . 

□  Analysis  Settings/Tool  name 

■  CFDTool:  The  CFDTool  is  the  default  fluidics  solver.  It  has  been  upgraded  in  MEMCAD  4.6  to  work 
with  a  Y2K-compatible  solver. 

■  ViewResult:  Bypass  the  setup  screens  and  view  previous  results  of  a  solution.  The  path  to  the  solution 
directory  is  set  from  the  Model  window. 

□  Analysis  Settings/Module 

■  NetFlow-C:  This  module  solves  the  containment  problem.  It  allows  only  for  the  solution  of  an 
advection-diffusion  equation.  Transport  is  by  diffusion  only. 

■  NetFlow-T:  This  module  solves  the  transport  problem.  It  enables  the  solution  of  the  full  Navier-Stokes 
equations  with  the  appropriate  flux  terms  for  the  various  transport. 

□  Analysis  Settings/Species 

Specifies  the  number  of  additional  species  present  in  the  flow,  including  the  carrier  fluid.  Choices  are  Carrier, 
Carrier+1,  and  Carrier+2  (0, 1,  and  2  species  present,  respectively). 

□  Analysis  Settings/Dimension 

■  3D:  Tells  the  solver  to  run  a  three  dimensional  solution  with  dimensions  X,Y,  and  Z. 

■  2D-Axi:  Allows  the  user  to  run  a  two  dimensional  axi-symmetric  solution.  The  model  has  to  be 
generated  in  the  XY-plane  and  should  be  1  element  thick  in  the  Z  direction.  The  X-axis  is  the  axis  of 
symmetry  and  the  Y-axis  is  the  radial  coordinate. 

■  2D:  allows  the  user  to  run  a  two  dimensional  solution.  The  model  has  to  be  generated  in  the  XY-plane 
and  should  be  1  element  thick  in  the  Z  direction. 

□  Analysis  Settings/ Analysis 

Sets  the  type  of  analysis  to  be  performed,  depending  on  the  concentration  of  the  fluid.  Dilute  and  Non-Dilute 
are  the  choices. 

□  OK 

Close  Tool  window;  return  control  to  MEMCAD  4. 
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NetFlow  -  Model Jnfo  hierarchical  window 

This  Model_Info  window  allows  access  to  the  names,  volumeMaterials,  patches,  and  species  windows  for  setting  up 
the  parameters  needed  by  the  solver. 


Each  button  opens 
a  window 

(described  separately). 


□  names  button 

Open  a  names  window  (see  page  F2-113). 

□  volumeMaterials  button 

Open  a  volumeMaterials  window  (see  page  F2-114). 

□  patches  button 

Open  a  patches  window  (see  page  F2-114). 

□  species  button 

Opens  a  species  window  (see  page  F2-115).  Note  that  this  button  does  not  appear  if  the  Tool  window  is  set  for 
a  Species=Carrier  value. 

□  Continue 

Close  the  window  and  open  the  boundary  conditions  menu  window  (see  page  F2-118). 

□  Cancel 

Close  window  and  return  to  MEMCAD  4  main  menu. 
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NetFlow  -  Model  Info  /  names  window 


The  names  window  allows  user  to  specify  the  Model  and  Material  Types.  This  information  is  placed  in  the  mbif  file 
when  it  is  translated. 
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□  Model  Type 

The  names  window  contains  a  line  for  each  FEM  created  during  the  MemBuilder  step.  User  can  specify  each 
FEM  to  be  either  movable,  fixed,  or  fluid. 

■  Movable:  If  set  to  movable,  the  NetFlow  solver  treats  the  part  as  a  solid;  this  means  that  the  thermal 
properties  of  the  solid  are  taken  into  account. 

■  Fluid:  Parts  assigned  a  fluid  type  have  flow  velocities,  pressures,  and  temperatures  as  degrees  of 
freedom.  Only  one  fluid  type  can  be  used  (such  as  water  or  air,  but  not  both).  Note  that  FlumeCAD 
includes  special  solvers  (BubbleSim  and  DropSim,  both  documented  in  this  volume),  which  are 
dedicated  to  solving  the  gas-liquid-solid  phase.  Such  problems  involve  surface  tension  and  can  be  used 
to  study  filling  of  structures,  bubble  transport,  and  fluid  ejection. 

■  Fixed:  A  fixed  part  type  will  be  not  enter  in  to  the  NetFlow  solver  calculations,  and  the  thermal 
properties  of  the  part  will  not  be  taken  into  account. 

□  Material  Type 

The  Material  Type  setting  is  not  used  for  the  NetFlow  solver.  See  the  User  Handbook  (page  U2-50)  to  learn 
how  to  use  dielectric  material  types  during  Convert  Model  to  obtain  patches  on  the  boundaries  between  solids 
and  fluids. 

□  Material  Info 

The  display-only  column  lists  the  current  material  assignment  for  the  FEM  line.  The  column  may  display  the 
material  assigned  during  3D  model  creation,  or  reassigned  from  the  Material  Property  Database,  depending  on 
when  the  window  is  viewed. 

□  OK 

Close  window  and  return  to  hierarchical  window  control. 
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NetFlow  -  Model  Jnfo  /  volumeMaterials  window 

This  volumeMaterials  window  allows  user  to  set  the  mechanical  constants  and  parameters  for  the  material  being  mod¬ 
eled  for  NetFlow.  The  window  is  the  same  as  the  DeviceCreation  /  Model  Setup  volumeMaterials  window  (see  the  Ref¬ 
erence  Manual,  page  R2-45). 
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NetFlow  -  ModeIJnfo  /  patches  window 

The  patches  window  allows  user  to  specify  names  for  the  patches  used  in  the  model.  The  window  is  the  same  as  the 
Device  Creation  /  Model  Setup  patches  window  (see  the  Reference  Manual,  page  R2-46). 
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NetFlow  -  Modeljnfo  /  species  window 

The  Species  window  allows  user  to  specify  the  species  characteristics  used  for  the  simulation.  The  species  present  are 
in  infinitely  dilute  form,  therefore  material  properties  are  not  needed.  The  Tool  window  must  have  Species  set  to  Car¬ 
rier +1  or  Carrier +2  to  access  this  window. 


□  Species 

This  view-only  window  lists  the  species  parameters  that  user  can  specify  for  simulation: 

■  Molecular  Weight:  scalar  measured  in  atomic  mass  units. 

■  Diffusion:  diffusion  constant  of  the  fluid,  measured  in  gm2/s. 

■  Mobility:  mobility  of  the  fluid,  measured  in  |im2/V*s. 

■  Conductivity:  conductivity  of  the  fluid,  measured  in  pS/gm. 

□  Species_l,  Species_2 

Designate  the  type  of  parameter  to  be  specified.  An  accompanying  Edit  window  (see  page  F2-116)  stores  the 
parameter  fields. Choices  include: 

■  Constant-Scalar:  Constant  with  a  single  value  over  the  entire  field. 

■  Polynomial-T:  Polynomial  function  of  temperature.  Six  constant  coefficients  are  allowed  to  enable  the 
definition  of  up  to  a  fifth  order  polynomial. 

■  Polynomial-T, VO:  Polynomial  function  of  temperature,  with  the  V0  term  used  to  define  an  initial  start¬ 
up  voltage  below  which  electrophoretic  effects  are  absent. 

■  Table-T:  Tabular  form  for  defining  temperature-property  pairs.  Linear  interpolation  is  used  for 
intermediate  values. 

□  Vary_l,  Vary_2 

Allows  specification  of  a  variable  tag  for  Simulation  Manager  runs. 
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NetFlow  -  Modeljnfo  /  species  settings  /  Edit  windows 

Each  LoadValue  type  selected  from  the  Species  window  uses  its  own  format  for  entering  parameters.  The  windows 
are  described  in  this  section. 

□  Constant-Scalar 

Constant  with  a  single  value  over  the  entire  field. 
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□  Polynomial-T 

Polynomial  function  of  temperature.  Six  constant  coefficients  are  allowed  to  enable  the  definition  of  up  to  a 
fifth  order  polynomial. 
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□  Polynomial-T,  VO 

Polynomial  function  of  temperature,  with  the  V0  term  used  to  define  an  initial  start-up  voltage  below  which 
electrophoretic  effects  are  absent. 
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□  Table-T 

Tabular  form  for  defining  temperature-property  pairs.  Linear  interpolation  is  used  for  intermediate  values. 
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NetFlow  -  hierarchical  Boundary  Conditions  menu  window 

The  NetFlow  BCs  hierarchical  boundary  conditions  menu  window  allows  user  to  choose  windows  to  specify  load 
conditions  for  the  NetFlow  fluidic  solver. 


from  ModeIJnfo 
hierarchical  window 
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Cancel  |  Simulate  if 

Each  button  opens  a  window 
(described  separately).  Each 
opened  window  returns  control 
to  this  hierarchical  window. 


Cancel;  return  to  MEMCAD  4  main  menu  Start  NetFlow  solver 

Open  results  and  Visualizer  module  windows 


□  surfaceBCs 

Open  window  (see  page  F2-119)  to  set  surface  boundary  conditions,  such  as  velocity  and  pressure. 


□  volumeBCs 

Open  window  (see  page  F2-125)  to  set  boundary  conditions  that  affect  an  entire  volume,  such  as  temperature 
conditions. 


□  transients 

Open  window  (see  page  F2-128)  to  set  conditions  for  transient  calculations. 

□  speciesSurfBCs 

Open  window  (see  page  F2-128)  to  set  conditions  for  species  surface  boundary  conditions. 

□  speciesVoIBCs 

Open  window  (see  page  F2-130)  to  set  conditions  for  species  volume  boundary  conditions. 

□  Simulate 

Launch  NetFlow  solver. 


□  Cancel 

Do  not  proceed  with  NetFlow  solver;  return  to  MEMCAD  4  main  menu. 
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NetFlow  -  NetFlow  BCs  /  surface  boundary  conditions  window 


The  surfaceBCs  window  sets  patch  surface  boundary  conditions.  Fix  types  and  load  conditions  for  these  surfaces  can 
be  specified  in  eight  Node/Face  sets  within  this  window.  When  setting  up  conditions,  the  LoadValue  button  needs  to 
correspond  to  the  appropriate  boundary  condition.  Not  all  combinations  of  boundary  conditions  and  LoadValue  are 
possible. 


□  surfaceBCs 

This  display-only  column  allows  user  to  define  up  to  eight  boundary  condition  sets  for  the  model. 


□  Fix  Type 

Enables  pull-down  menu  with  many  load  type  choices: 

■  Wall:  This  represents  a  “no  slip”  boundary  condition  at  the  selected  patches.  The  no-slip  boundary 
condition  is  an  assumption  that  is  made  in  Newtonian  Fluids.  Newtonian  fluids  also  form  a  parabolic 
velocity  profile  in  steady  state.  No  LoadValue  is  applied. 

■  Velocity:  Possible  LoadValues  are  Velocity  and  Vector. 

■  X- Velocity:  This  option  specifies  a  uniform  velocity  in  the  X-direction.  The  option  is  very  useful  when 
applying  symmetric  boundary  conditions  for  which  the  velocity  perpendicular  to  the  axis  of  symmetry 
is  0.  A  Scalar  LoadValue  is  used. 
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■  Y- Velocity:  This  option  specifies  a  uniform  velocity  in  the  Y-direction.  The  option  is  very  useful  when 
applying  symmetric  boundary  conditions  for  which  the  velocity  perpendicular  to  the  axis  of  symmetry 
is  0.  A  Scalar  LoadValue  is  used. 

■  Z- Velocity:  This  option  specifies  a  uniform  velocity  in  the  Z-direction.  The  option  is  very  useful  when 
applying  symmetric  boundary  conditions  for  which  the  velocity  perpendicular  to  the  axis  of  symmetry 
is  0.  A  Scalar  LoadValue  is  used. 

■  Pressure:  A  constant  pressure  on  the  selected  patch  is  specified.  A  Scalar  LoadValue  is  used. 

■  Symmetry:  The  selected  patch  is  a  plane  of  symmetry.  It  is  important  that  the  entire  problem,  not  just 
the  geometry,  is  symmetric.  No  LoadValue  is  applied. 

■  Mass  Flux:  This  option  is  valid  only  for  the  second  and  third  species;  a  value  cannot  be  provided  for  the 
carrier  fluid.  A  Vector  LoadValue  is  used. 

■  FlowRate:  The  flow  rate  follows  the  sign  convention  of  all  MEMCAD  4  boundary  conditions.  A 
positive  flow  rate  is  along  the  normal  of  the  selected  patch.  Since  in  MEMCAD  4  the  normal  always 
points  to  the  volume,  the  FlowRate  is  into  the  patch  when  a  positive  value  is  specified.  A  Scalar 
LoadValue  is  used. 

■  Temperature:  Applies  a  constant  temperature  on  the  selected  patch.  A  Scalar  LoadValue  is  used. 

■  Retardation:  This  boundary  condition  allows  application  of  wall  effects  such  as  chemical  forces 
experienced  by  the  species  near  a  wall,  slowing  down  the  diffusion  and  changing  the  species  mobility 
near  the  wall.  A  Retardation  LoadValue  is  used. 

■  Voltage:  This  boundary  condition  applies  a  constant  potential  on  the  specified  patch.  Note  that  a 
potential  difference  must  be  specified;  therefore,  two  Voltage  BC  specifications  are  required.  A  Scalar 
LoadValue  is  used. 

■  Current:  This  boundary  condition  applies  a  current  through  the  specified  patch.  A  Scalar  LoadValue  is 
used. 

■  Wall  Effects:  This  boundary  condition  allows  the  specification  of  electroosmotic  effects  in  the  flow. 
EDL,  EO  _M  ability,  and  Partial_EO  LoadValues  can  be  used.  Slip  wall  effects  also  can  be  taken  into 
account  with  the  Slip  LoadValue. 

□  Patch  1,2,3 

Choose  from  any  of  the  patch  names  previously  set  up. 

□  andl,2,3 

Boolean  or  that  defines  intersection.  Use  the  Boolean  or  for  the  Fix  Type  to  apply  independently  to  each  of  the 
selected  patches.  The  Boolean  and  is  disabled  in  NetFlow. 
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□  Load  Value 

Choose  from  several  load  types.  Each  type  uses  a  different  edit  window;  the  entire  set  of  edit  windows  is 
shown  as  part  of  the  LoadValue  explanation. 

■  Scalar:  one-dimensional  quantity,  such  as  pressure  load  or  temperature. 


none  jc.5 
PotT)}  Cancel . 


■  Vector:  Used  with  the  Velocity  FixType.  This  defines  a  uniform  velocity  profile  in  the  direction  of  the 
vector  specified  on  the  selected  patch. 
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■  Velocity:  Used  with  the  Velocity  FixType.  This  is  a  quadratic,  linear  or  constant  polynomial  in  all  three 
directions  for  the  velocity  on  the  selected  patch.  For  each  velocity  Ux,  Uy  and  Uz  a  separate 
polynomial  can  be  given.  For  example,  Pouiselle  Flow  in  a  pipe  yields  the  following  profile: 


where  C  is  a  constant  and  Yq  is  the  pipe  radius. 

For  this  example,  the  Edit  Velocity  window  below  requires  that  only  the  first  U x  line  needs  to  be  filled 
in.  The  following  column  values  are  required  for  this  window: 

2 

ColumnA  =  C,  ColumnC  =  -C,  Column  =  1.  All  remaining  columns  are  zero. 

The  example  shown  is  filled  out  for  the  constant  C-5. 
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■  Retardation:  The  following  fields  are  provided  as  input. 

Max_distance:  The  distance  over  which  the  wall  effects  are  felt.  This  can  be  seen  as  a 
Wall-Effect-Boundary  layer. 

Diffusive:  Specifies  the  diffusive  retardation  coefficient;  i.e.,  a  mechanism  for  slowing  the  diffu¬ 
sion  of  the  species  in  the  region  within  the  Max_distance  parameter. 

Convective:  Specifies  the  convective  retardation  coefficient;  i.e.,  the  change  in  the  mobility  of  the 
species  within  the  near-wall  region. 

In  both  coefficient  cases  the  retardation  is  applied  as  a  smooth  function  from  a  maximum  near  the 
wall  to  a  zero  retardation  value  at  the  specified  Maxjiistance  value. 


Edit  Retardation 
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Cancel 
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■  Slip:  The  Mean  Free  Path  and  the  MAC  (Momentum  Accommodation  Coefficient)  are  used  in  the  slip 
equation  to  study  wall  effects  in  microchannels.  It  is  used  exclusively  with  the  Wall  Effects  FixType. 
The  slip  velocity  at  the  wall  can  be  expressed  as  follows: 

SlipVelocity  =  ok(  1 

V  / wall 

where  U  is  the  velocity,  n  is  the  direction  normal  to  the  velocity,  O  is  the  streamwise  momentum  accom¬ 
modation,  expressed  as  follows: 


where  CTm  is  the  momentum  accommodation  coefficient  (MAC).  Usually,  <Jm  =1  but  it  can  vary  within 
the  range  0  <  Om  <  1 . 


K  is  the  Knudsen  number,  which  can  be  expressed  as: 


K_  - 


MeanFreePath 


n  CharacteristicLength 
If  this  value  is  smaller  than  0.01  the  slip  effect  is  neglected. 
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■  EDL:  EDL  stands  for  Electrical  Double  Layer,  which  is  governed  by  the  Poisson-Boltzmann  equation. 
Four  values  are  needed:  Zeta-Potential,  Charge,  Ionic, Concentration,  and  Dielectric  coefficient  of  the 
carrier  fluid.  Setting  the  EDL  LoadValue  activates  the  solution  of  the  Poisson-Boltzmann  equation  for 
the  double  layer  to  couple  in  electroosmotic  effects.  This  setting  should  be  used  with  care,  since  the 
resolution  of  the  double  layer  in  most  conventional  problems  requires  a  very  dense  mesh  near  the  walls. 


■  EO_Mobility:  This  setting  enables  the  representation  of  the  electroosmotic  effects  on  the  walls  through 
a  mobility  that  defines  the  velocity  achieved  by  the  charged  carrier  at  the  edge  of  the  double  layer  in  the 
electric  field.  Note  that  this  implies  a  relaxation  of  the  no-slip  wall  boundary  condition,  since  the  near- 
wall  region  is  no  longer  solved.  The  EO_Mobility  is  a  positive  number  for  a  negatively  charged  wall. 
(In  other  words,  a  positive  value  will  result  in  flow  towards  the  negative  electrode.) 


■  Partial_EO:  This  setting  measures  the  EO  mobility  when  the  surface  properties  change  across  the 
length  of  a  channel.  To  use  this  LoadValue  with  a  specific  patch,  define  a  volume  where  the  surface 
changes  occur  with  the  first  six  entries,  and  enter  the  new  EO  mobility  for  this  region  with  the  last  entry. 
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□  Variable 

Setting  for  Simulation  Manager.  A  Fixed  setting  (the  default)  uses  constant  Load  Values,  which  are  not 
changed  during  a  Simulation  Manager  run.  When  a  Load  Value  is  to  be  a  variable  modified  by  a  Simulation 
Manager  trajectory,  one  of  eight  NF_BC  variables  can  be  used,  corresponding  to  the  lines  in  the  netFlowBC 
Simulation  Manager  setup  window. 

□  Transient 

Setting  for  transients  surface  boundary  condition.  A  Fixed  setting  (the  default)  does  not  apply  transients. 
When  a  Node  Set  is  to  be  modified  by  a  transient  event,  one  of  two  Transient  variables  are  available  for 
assignment.  The  variable  is  defined  in  the  transients  window  (see  page  F2-128).  With  this  technique  a  peri¬ 
odic  pressure  or  input  flow  can  be  generated.  By  setting  the  Transient  variable  the  input  on  that  node  set  is 
coupled  to  the  waves  or  function  defined  in  the  transients  window. 
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NetFlow  -  NetFlow  BCs  /  volume  boundary  conditions  window 

The  volumeBCs  window  sets  volume  boundary  conditions.  Types  and  load  conditions  for  these  volumes  can  be  spec¬ 
ified  in  Sets  within  this  window. 


□  volumeBCs 

This  display-only  column  allows  user  to  define  up  to  eight  boundary  condition  sets  for  the  model. 

□  BCiype 

Enables  pull-down  menu  with  several  choices: 

■  none:  Do  not  apply  a  boundary  condition  to  the  volume. 

■  Mass  Flow:  Specify  a  mass  flow  in  three  directions — Mx,  My,  Mz.  The  Mass  Flow  LoadValue  is  used. 

■  Velocity:  Apply  a  velocity  to  the  selected  part.  A  Vector,  Velocity,  or  InputFile  LoadValue  is  used. 

■  Heat  Generation:  Heat  can  be  generated  in  a  part.  The  part  is  usually  a  solid  and  not  a  fluid.  A  Scalar 
LoadValue  is  used. 

■  Temperature:  Apply  a  fixed  temperature  to  the  selected  part.  A  Scalar  LoadValue  is  used. 

■  E-Field:  The  E-Field  in  the  volume  can  be  read  from  a  pre-computed  mbif  file.  Use  the  InputFile 
LoadValue  and  specify  the  directory  path. 

■  Temperature  Cycle:  Apply  a  temperature  cycle  to  the  part. 

■  Joule  Heat:  Includes  Joule  heating  effects.  Voltage/Current  BCs  must  be  set  and  the  carrier  should  have 
an  electrical  conductivity  material  property  defined. 

■  T-Inital:  specifies  initial  temperature  conditions. 
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□  Part 


Enables  pull-down  menu  with  a  list  of  all  the  parts  modeled  and  stored  in  the  mbif  file. 

□  Load  Value 

Choose  from  several  load  types.  Each  type  uses  a  different  edit  window,  shown  below  with  the  LoadValue 
description. 

■  Scalar:  Used  to  specify  a  one-dimensional  value,  such  as  Temperature  or  Heat  Generation. 


■  Mass  Flow:  Used  to  specify  a  three-dimensional  value  for  Mass  Flow. 


■  Velocity:  Used  with  the  Velocity  BCType.  This  is  a  quadratic,  linear  or  constant  polynomial  in  all  three 
directions  for  the  velocity  on  the  selected  part.  For  each  velocity  Ux,  Uy,  and  Uz  a  separate  polynomial 
can  be  given. 


■  PreCompute:  Decouples  the  momentum  and  species  equations.  The  velocity  is  pre-computed  and  an 
advection-diffusion  analysis  for  the  species  follows.  No  LoadValue  edit  window  is  used. 

■  Pre-Compute  File:  Specify  a  result  mbif  file  for  the  pre-compute  values.  See  Input  File  for  window 
description. 
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■  Input  File:  Specify  a  result  mbif  file  from  another  simulation  as  the  source  data  input  initial  conditions. 
The  default  directory  is  the  directory  from  which  MEMCAD  is  launched. 
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■  Joule_Transient:  Specify  a  scalar  value  for  a  joule  transient. 
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□  Variable 

Setting  for  Simulation  Manager.  A  Fixed  setting  (the  default)  uses  constant  Load  Values,  which  are  not 
changed  during  a  Simulation  Manager  run.  When  a  Load  Value  is  to  be  a  variable  modified  by  a  Simulation 
Manager  trajectory,  one  of  eight  NF_BC  variables  can  be  used,  corresponding  to  the  lines  in  the  netFlowBC 
Simulation  Manager  setup  window. 

□  Transient 

Setting  for  transients  surface  boundary  condition.  A  Fixed  setting  (the  default)  does  not  apply  transients. 
When  a  Node  Set  is  to  be  modified  by  a  transient  event,  one  of  two  Transient  variables  are  available  for 
assignment.  The  variable  is  defined  in  the  transients  window  (see  page  F2-128). 
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NetFlow  -  NetFlow  BCs  /  species  surface  boundary  conditions  window 

The  speciesSurfBCs  window  sets  surface  boundary  conditions  for  the  species.  Fix  types  and  load  conditions  for  these 
surfaces  can  be  specified  in  eight  Node/Face  sets  within  this  window.  When  setting  up  conditions,  the  LoadValue 
button  needs  to  correspond  to  the  appropriate  boundary  condition.  Not  all  combinations  of  boundary  conditions  and 
LoadValues  are  possible. 


□  speciesSurfBCs 

This  display-only  column  allows  user  to  define  up  to  eight  boundary  condition  sets  for  the  model. 

□  Species 

Can  choose  either  Species  1  or  Species  2  for  the  setting. 

□  Patch  1,2, 3 

Choose  from  any  of  the  patch  names  previously  set  up. 

□  andl,2,3 

Boolean  and/or  that  defines  union  or  intersection.  Use  the  Boolean  or  for  the  Fix  Type  to  apply  independently 
to  each  of  the  selected  patches.  Use  the  Boolean  and  for  the  Fix  Type  to  apply  to  the  common  intersection  of 
the  selected  patches. 

□  Load  Value 

Choose  from  several  load  types.  Each  type  uses  a  different  edit  window;  the  entire  set  of  edit  windows  is 
shown  as  part  of  the  LoadValue  explanation. 

■  Scalar:  a  fixed  value  on  a  patch. 
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■  Inj_Gaussian:  the  injection  of  a  finite  width  Gaussian  plug  through  the  patch.  The  Full- Width  at  Half 
Max  (FWHM)  value  can  be  defined  in  seconds  or  in  microns. 
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■  Inj_Plug:  the  injection  of  a  finite  width  square  plug  through  the  patch.  The  Length  of  the  plug  can  be 
defined  in  seconds  or  in  microns. 
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□  Variable 

Setting  for  Simulation  Manager.  A  Fixed  setting  (the  default)  uses  constant  Load  Values,  which  are  not 
changed  during  a  Simulation  Manager  run.  When  a  Load  Value  is  to  be  a  variable  modified  by  a  Simulation 
Manager  trajectory,  one  of  eight  NF_BC  variables  can  be  used,  corresponding  to  the  lines  in  the  netFlowBC 
Simulation  Manager  setup  window. 

□  Transient 

Setting  for  transients  surface  boundary  condition.  A  Fixed  setting  (the  default)  does  not  apply  transients. 
When  a  Node  Set  is  to  be  modified  by  a  transient  event,  one  of  two  Transient  variables  are  available  for 
assignment.  The  variable  is  defined  in  the  transients  window  (see  page  F2-128).  With  this  technique  a  peri¬ 
odic  pressure  or  input  flow  can  be  generated.  By  setting  the  Transient  variable  the  input  on  that  node  set  is 
coupled  to  the  waves  or  function  defined  in  the  transients  window. 
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NetFlow  -  NetFlow  BCs  /  species  volume  boundary  conditions  window 

The  species VolBCs  window  sets  volume  boundary  conditions  for  the  species.  Types  and  load  conditions  for  these 
volumes  can  be  specified  in  Part  sets  within  this  window. 
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□  speciesVolBCs 

This  display-only  column  allows  user  to  define  up  to  eight  boundary  condition  sets  for  the  model. 

□  Species 

Can  choose  either  Species  1  or  Species  2  for  the  setting. 

□  BCType 

Enables  pull-down  menu  with  several  choices: 

■  none:  does  not  apply  a  boundary  condition  to  the  volume. 

■  Fixed:  fixes  the  selected  volume. 

■  Initial:  This  boundary  condition  is  used  only  for  transient  flow  analysis.  A  patch  needs  to  be  associated 
and  selected  in  order  to  pose  a  numerically  correct  problem.  The  selected  method  enters  Mass  Fraction 
into  the  flow  as  desired.  The  transient  problem  needs  a  Dirichlet  BC  for  the  concentration,  which  is 
chosen  to  be  0  on  the  selected  patch.  Thus,  the  selected  patch  should  be  far  enough  upstream  in  order  for 
the  Mass  Fraction  not  to  reach  the  patch  at  any  time. 

□  Part 

Enables  pull-down  menu  with  a  list  of  all  the  parts  modeled  and  stored  in  the  mbif  file. 

□  Load  Value 

Choose  from  several  load  types.  Each  type  uses  a  different  edit  window,  shown  below  with  the  LoadValue 

description. 

■  Scalar:  Used  to  specify  a  one-dimensional  value,  such  as  Temperature  or  Heat  Generation. 
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■  Location:  Confines  the  concentration,  given  by  the  Mass  Fraction,  to  the  selected  minimum  and 
maximum  values  of  X,  Y,  and  Z.  MEMCAD  4  finds  all  nodes  that  are  within  these  bounds  and  applies 
the  Mass  Fraction.  If  the  minimum  and  maximum  in  any  direction  are  left  unchanged  (i.e  zero),  it  means 
that  “no  filter”  (all  elements  between  -  °°  and  +  °°)  is  applied  in  this  direction.  Using  this  LoadValue  can 
lead  to  solver  uncertainties,  due  to  the  sharp  edges  created  by  the  Location  setting.  The  Gaussian  setting 
is  preferred. 
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Gaussian:  Applies  a  Gaussian  distribution  according  to  the  following  formula: 

A  —  Aq 


A  exp 


Gx 
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J  LV 


where  A  is  the  amplitude,  and  0  is  the  change  in  a  particular  direction.  Typically,  only  one  dimension  is 
given  a  non-zero  0  value,  which  results  in  a  plug  of  species  given  a  Gaussian  distribution  in  one  axial 
direction.  The  Gaussian  setting  greatly  reduces  solver  instability,  and  is  the  preferred  LoadValue  for 
most  calculations. 
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Input  File:  Specify  a  result  mbif  file  from  another  simulation  as  the  source  data  input  boundary 
conditions.  The  default  directory  is  the  directory  from  which  MEMCAD  is  launched. 
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□  Variable 

Setting  for  Simulation  Manager.  A  Fixed  setting  (the  default)  uses  constant  Load  Values,  which  are  not 
changed  during  a  Simulation  Manager  run.  When  a  Load  Value  is  to  be  a  variable  modified  by  a  Simulation 
Manager  trajectory,  one  of  eight  NF_BC  variables  can  be  used,  corresponding  to  the  lines  in  the  netFlowBC 
Simulation  Manager  setup  window. 

□  Transient 

Setting  for  transients  surface  boundary  condition.  A  Fixed  setting  (the  default)  does  not  apply  transients. 
When  a  Node  Set  is  to  be  modified  by  a  transient  event,  one  of  two  Transient  variables  are  available  for 
assignment.  The  variable  is  defined  in  the  transients  window  (see  page  F2-134). 
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NetFlow  -  NetFlow  BCs  /  transients  boundary  conditions  window 

The  transients  window  enables  the  user  to  specify  boundary  conditions  associated  with  transient  calculations.  Tran¬ 
sient  variables  can  be  assigned  in  the  surfaceBC  or  volumeBC  window,  and  the  type  of  transient  response  to  use  for 
the  problem  can  be  specified  in  this  window.  Two  different  transient  types  can  be  specified.  For  any  transient  analy¬ 
sis,  the  appropriate  parameters  should  be  set  up  in  the  NetFlow  Tool  window. 


Transients  are  restricted  to  the  following  boundary  condition  settings:  SurfaceBCs — Velocity,  X,Y,Z 
Velocity,  Pressure,  FlowRate  Temperature,  Heat  Flux,  Concentration',  VolumeBCs — Heat  Generation. 
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□  transients 

This  display-only  column  allow  user  to  define  up  to  eight  boundary  condition  transient  sets  for  the  model. 

□  ITansient 

This  column  allows  assignment  of  one  of  two  transient  curve  variables:  Transient  1  or  Transient 2.  The  equiva¬ 
lent  TransientBCl  and  TransientBC2  variables  are  assigned  in  either  the  surfaceBC  or  volumeBC  window  to  a 
desired  parameter  for  transient  analysis. 

□  Curve  Type  /  Edit 

The  NetFlow  transient  capability  allows  one  of  four  different  transient  waveforms  to  be  specified  for  the  anal¬ 
ysis.  The  specific  waveform  characteristics  are  specified  in  the  Edit  window  accessed  from  the  Edit  button. 
These  Edit  windows  are  further  described  on  page  F2-134. 

■  Square:  Assigns  a  square  wave  as  the  transient  curve.  The  square  wave  can  be  controlled  by  adjusting 
its  period,  split  (duty  cycle),  or  amplitude. 

■  Periodic:  Assigns  a  periodic  curve  as  the  transient.  The  periodic  curve  can  be  controlled  by  adjusting  its 
period,  mean,  or  amplitude. 

■  Curve:  Assigns  a  table  of  time  step  and  property  values  in  order  to  create  a  customized  curve. 

■  Sawtooth:  Produces  a  sawtooth  waveform. 
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□  Variable 


Setting  for  Simulation  Manager.  A  Fixed  setting  (the  default)  uses  constant  Load  Values,  which  are  not 
changed  during  a  Simulation  Manager  run.  When  a  Load  Value  is  to  be  a  variable  modified  by  a  Simulation 
Manager  trajectory,  one  of  two  FluidBC  variables  can  be  used,  corresponding  to  the  lines  in  the  fluidBCs  Sim¬ 
ulation  Manager  setup  window. 


NetFlow  -  NetFlow  BCs  /  transients  BC  window  /  Edit  CurveType  windows 

Each  transient  LoadValue  type  selected  from  the  Edit  CurveType  window  uses  its  own  format  for  entering  parame¬ 
ters.  The  windows  are  described  in  more  detail  in  this  section. 

□  Square 

The  square  wave  can  be  controlled  by  adjusting  its  period,  split  (duty  cycle),  or  amplitude.  Period  is  the  time 
between  the  start  of  two  square  waves.  Split  is  a  value  between  0  and  1 .  If  Split  is  0.5,  the  square  wave  is  at  its 
maximum  for  half  a  period.  If  the  Split  is  0.3,  then  the  square  wave  is  at  its  maximum  for  30%  of  the  period. 
Amplitude  is  the  maximum  value  of  the  square  wave. 
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□  Periodic 

The  periodic  wave  can  be  controlled  by  adjusting  its  period,  mean,  or  amplitude:  Period  is  the  time  between 
the  start  of  two  periodic  waves.  Mean  moves  the  periodic  curve  vertically  (up  or  down).  Amplitude  is  the  max¬ 
imum  value  of  the  periodic  wave. 
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□  Curve 

The  custom  curve  is  described  with  a  table  of  time  step  and  property  values.  Users  can  enter  up  to  10  different 
sets  of  time  and  amplitude  values  to  simulate  almost  any  type  of  transient  curve. 


□  Sawtooth 

The  sawtooth  is  described  with  a  combination  of  time  and  temperature  values,  along  with  a  slope  specification 
for  the  sawtooth  angle. 
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NetFlow  -  hierarchical  solutions  menu  window 

The  fluidic  results  hierarchical  solutions  menu  window  allows  user  to  choose  windows  to  view  result  tables  and  mbif 
file  paths  used  during  the  computation. 


appears  when  NetFlow 
solution  is  complete 
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—  Dens  ! 

Return  to  MEMCAD  4  main  menu 


Each  button  opens  a  window 
(described  separately).  Each 
opened  window  returns  control 
to  this  hierarchical  window. 


□  fluidDomain 

Open  window  (see  page  F2-137)  to  display  table  of  fluid  results. 

□  spIntegral/Table 

Open  window  (see  page  F2-138)  to  display  table  of  splntegral  results. 

□  spIntegral/Graph 

Open  window  (see  page  F2-138)  to  display  graph  of  splntegral  results. 

□  file  Array 

Open  window  (see  page  F2-139)  to  display  directory  paths  of  all  saved  mbif  Files. 

□  files 

Open  window  (see  page  F2-139)  to  display  directory  path  of  solution  results. 

□  Done 

Close  window  and  return  control  to  MEMCAD  4. 

□  File/Print 

Enable  printing  of  the  individual  window  results  to  a  single  file.  The  text  format  file  can  be  used  to  interface  to 
other  programs.  The  file  is  located  in  the  directory  from  which  MEMCAD  4  starts. 
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NetFlow  -  Results  /  fluid  solution  window 

This  fluidDomain  window  displays  the  results  from  the  NetFlow  solver. 
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□  fluidDomain 

List  the  type  of  velocity  for  the  solver  solution.  Maximum  velocity  as  well  as  the  individual  X,Y,  and  Z  com¬ 
ponents  are  displayed. 

□  Maximum 

The  maximum  velocity  in  the  channel  calculated  by  the  solver.  The  value  shown  is  an  absolute  value. 

□  Minimum 

The  minimum  velocity  in  the  channel  calculated  by  the  solver.  The  sign  of  the  value  is  derived  from  the  nor¬ 
mal  vector;  a  minimum  value  can  represent  maximum  velocity.  If  the  value  in  this  column  represents  the  max¬ 
imum  velocity  for  the  model,  the  same  value  will  appear  in  the  Maximum  column. 

□  OK 

Close  window. 
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NetFlow  -  Results  /  spintegral  Table  and  Graph  solution  windows 

These  windows  display  the  numerical  and  graphical  species  mass  results  from  the  simulation.  The  transient  time  and 
output  step  settings  for  the  table  and  graph  are  set  in  the  NetFlow  Tool  window  (see  page  F2-109). 
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□  spintegral  Time 

This  column  and  abscissa  in  the  graph  display  the  output  step  intervals  set  in  the  NetFlow  Tool  window.  In  the 
above  example,  output  mbif  Fdes  are  written  every  0.2  seconds,  and  the  mass  results  are  displayed. 

□  Mass 

This  column  and  ordinate  display  the  species  mass  in  the  carrier  fluid. 

□  OK 

Close  window. 
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NetFlow  -  Results  /  fileArray  window 

The  fileArray  window  lists  the  directory  paths  for  each  output  mbif  file  written  during  the  simulation.  The  number  of 
files  is  determined  by  the  output  time  step  length  and  total  transient  simulation  run  time,  as  set  in  the  SwitchSimTool 
window  (see  page  F2-109). 
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NetFlow  -  Results  /  result  files  window 

This  files  window  displays  the  path  to  the  result  mbif  file. 
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Simulation  Manager  controls  for  NetFlow 

When  NetFlow  boundary  conditions  variables  are  set  up,  the  Simulation  Manager  can  assign  a  trajectory  to  the  vari¬ 
able  and  produce  a  batch  run  simulation,  creating  separate  mbit  files  for  each  trajectory  step.  This  section  lists  the 
NetFlow-specific  controls  accessed  from  the  Simulation  Manager.  For  a  complete  reference  description  of  general 
Simulation  Manager  capability,  refer  to  the  Reference  Manual,  page  R2-94. 


□  View  NF_BC 


This  gray  button  provides  convenient  viewing  access  for  checking  the  NetFlow  boundary  condition  windows, 
which  are  described  in  detail  on  pages  page  F2-116  through  page  F2-135.  The  windows  are  read-only  and 
cannot  be  modified  from  this  access  point. 

□  netFlowBCs 

The  green  button  opens  a  variable  assignment  window,  which  allows  the  user  to  define  up  to  eight  fluid 
boundary  conditions.  The  variables  arc  assigned  in  the  species,  surfaceBC,  speciesSurfBCs,  volumeBC,  and 
species VolBCs  windows,  and  equated  in  this  netFlowBC  window  to  the  trajectory  type  set  for  the  Simulation 
Manager  problem.  The  functionality  of  the  window  is  similar  to  the  SimMan/MechBCs  window  described  in 
the  Reference  Manual,  page  R2-110. 

□  Simulate 

Starts  the  Simulation  Manager  calculations.  When  finished,  a  hierarchical  results  window  appears,  and  allows 
access  to  a  velocity  table  and  accompanying  graph. 
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Visualizer  -  Query  function 

Clicking  on  the  Query  button  on  the  Visualizer  Control  Panel  opens  a  set  of  hierarchical  windows.  Through  these 
windows,  the  user  can  select  surface  or  volume  extraction  criteria  for  result  values  stored  in  mbif  files.  An  executable 
then  performs  the  extraction  and  writes  an  output  file  for  display  or  input  into  another  solver. 
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_►  Each  button  opens  a  window 
(described  separately).  Control 
-►  returns  to  this  window. 


Begins  Query  simulation;  returns 
to  the  hierarchical  Simulation 
Results  window  when  finished. 


□  On  a  Patch 

Opens  window  that  allows  user  to  set  query  parameters  on  a  patch.  See  page  F2-142  for  further  details. 

□  Between  Patches 

Opens  window  that  allows  user  to  set  query  parameters  between  two  patches.  See  page  F2-143  for  further 
details. 

□  In  a  Volume 

Opens  window  that  allows  user  to  set  query  parameters  on  a  part  for  a  user-defined  volume. 

See  page  F2-144  for  further  details. 

□  In  a  Region 

Opens  window  that  allows  user  to  set  query  parameters  on  a  part  for  a  point,  volume,  or  an  ellipse  within  a 
user-defined  volume.  See  page  F2-145  for  further  details. 

□  Cancel 

Query  is  not  performed,  and  control  is  returned  to  the  Visualizer  Control  function  bar. 

□  Execute 

Launches  Query  solver. 
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Query  -  On  a  Patch 

The  patch  query  allows  users  to  extract  result  values  for  a  specified  patch  within  a  MEMCAD  model  mbif  file.  Not  all 
field  types  are  applicable;  the  relevant  types  are  listed  below. 


IIS 

W3E3I 

JfQump  On  A  Patch 

HhoT 

;  ",-"Ti  n 

- . 

Fils  Help 

s  wn  *  “ 

”  ■  quiiyjsbti 

Label 

:  Query_Label 

IpatchQueiy 

iPi.  :|  m 

Pitch  Qubif  | 

Patch  |  Field 

.  Query!  j 

none 

j^Jjnone 

d; 

vrfi  rut 

Query2 

none 

j?J)none 

d 

QueiyG  i 

none 

_ijjnore 

3 

in  nun  >1 

]  Query4  j 

none 

^ilnore 

d 

_ 

!  QueryS  1 

none 

^llnone 

— 3 

ten  iICjsiio 

•  Queiy6  j 

none 

^llnone 

iii 

j  Query7  j 

none 

-llnone 

di 

»:  wwti 

t»CUt€ 

;  Queiy9  j 

none 

jjjnone 

■ 

dm 

Cancel  | 

i 


1  none 

o\ 

HeatFluK 

Current_Fiux 

Charge 

— 

Force 

Species_l_FlLrx 

Species  JLFlux 

SpeciesJII_Flu* 

* 

□  Label 

An  editable  field  that  allows  user  to  change  the  label  name.  The  label  will  appear  in  the  results  window. 

□  Patch  Query 

Indicates  the  query  number.  Each  set  can  have  a  unique  patch  name  and/or  field  type  setting. 

□  Patch 

User  selects  a  patch  name  from  a  pull-down  menu. 

□  Field 

User  selects  a  calculation  type  from  a  pull-down  menu. 

■  Heat_Flux:  calculates  heat  flux  through  the  selected  patch. 

■  Current_Flux:  calculates  current  flux  through  the  selected  patch. 

■  Species_/V_Flux:  calculates  species  flux  through  a  selected  patch. 

■  AvgVoltage:  calculates  average  voltage  on  a  patch. 

■  Avg,  Max,  MinTemperature:  calculates  average,  maximum,  or  minimum  temperature  on  a  patch, 
respectively. 

■  FlowRate:  measures  the  flow  rate  through  a  patch. 
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Query  -  Between  Patches 

The  between  patches  query  allows  users  to  extract  result  values  dependent  on  two  patches  within  a  MEMCAD  model 
mbif  file. 


□  Label 

An  editable  field  that  allows  user  to  change  the  label  name.  The  label  will  appear  in  the  results  window. 


□  Relative  Query 

Indicates  the  query  number.  Each  set  can  have  unique  patch  names  and/or  field  type  setting. 

□  Patchl,  2 

User  selects  a  patch  name  from  a  pull-down  menu.  Both  patch  fields  must  be  selected  in  order  to  complete  the 
query. 

□  Field 

User  selects  a  calculation  type  from  a  pull-down  menu. 

■  Resistance:  calculates  resistance  between  the  two  selected  patches. 

■  Voltage:  calculates  voltage  between  the  two  selected  patches. 

■  PressureDrop:  measures  change  in  pressure  between  the  two  selected  patches. 
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Query  -  In  a  Volume 

The  volume  query  allows  the  user  to  extract  result  values  for  a  specified  volume  within  a  MEMCAD  model  mbif  file. 
Not  all  field  types  are  applicable;  the  relevant  types  are  listed  below. 
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□  Label 

An  editable  field  that  allows  user  to  change  the  label  name.  The  label  will  appear  in  the  results  window. 

□  Name  Query 

Indicates  the  query  number. 

□  Name 

User  selects  a  part  name  from  a  pull-down  menu.  A  field  calculation  type  is  applied  to  this  name. 

□  Field 

User  selects  a  calculation  type  from  a  pull-down  menu. 

■  Species/V_Integral:  determines  the  percentage  of  defined  species  that  passes  through  a  part. 

■  Max,  MinTemperature:  calculates  maximum  or  minimum  temperature  on  the  defined  part. 
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Query  -  In  a  Region 

The  region  query  allows  users  to  extract  result  values  for  a  portion  of  a  specified  volume  within  a  MEMCAD  model. 
The  region  can  be  a  point  of  an  inner  volume  within  a  specific  volume.  A  point  is  specified  with  a  single  XYZ  coor¬ 
dinate  defining  a  location.  An  inner  volume  is  specified  with  two  XYZ  coordinates  defining  the  volume  space. 


Clicking  on  the  Region  Query  button  opens  a  setup  window. 

□  Label 

An  editable  field  that  allows  user  to  change  the  label  name.  The  label  will  appear  in  the  results  window. 

□  Region  Query 

Indicates  the  query  number. 

□  Name 

User  selects  a  part  name  from  a  pull-down  menu. 
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□  Location 


User  selects  a  point,  volume,  or  elliptical  region  from  a  pull-down  menu  to  be  defined.  Each  selection  opens  a 
window  with  editable  values.  The  choices  are  Point ,  Volume ,  or  Ellipse. 

■  Point:  user  defines  X,  Y,  and  Z-coordinates  of  the  point. 


■  Volume:  user  defines  volume  space  and  a  mesh  for  query  calculation.  If  any  pair  of  min/max  values  is 
left  at  0,  no  filter  is  applied. 


XMin,  XMax:  defines  the  X-coordinates  for  the  volume. 

YMin,  YMax:  defines  the  Y-coordinates  for  the  volume. 

ZMin,  ZMax:  defines  the  Z-coordinates  for  the  volume. 

NX,  NY,  NZ:  Defines  number  of  mesh  elements  for  query  calculation  in  the  X,  Y,  and  Z-dircc- 
tions,  respectively. 


Query  mesh 

The  NX,  NY,  and  NZ  mesh  elements  defined  in  this  window  are  calculated  along  with  the  mesh  defined 
in  the  solver.  To  preserve  maximum  accuracy,  keep  the  query  mesh  density  as  close  as  possible  to  the 
solver  mesh  density. 
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■  Ellipse:  User  defines  boundaries  of  the  ellipse  and  a  mesh  for  query  calculation.  This  is  typically  used 
when  firing  a  laser  through  a  channel  to  determine  the  amount  of  species  present  in  a  specific  volume. 


xO,  yO,  zO:  defines  the  origin  of  the  ellipse. 

x_size,  y_size,  z_size:  the  value  assigned  in  each  direction  is  the  diameter  of  the  laser  where  the 
intensity  is  greater  than  or  equal  to  e"1. 

Nel_x,  Nel_y,  Nel_z:  Defines  the  mesh  in  the  x,  y,  and  z-directions.  Keep  this  mesh  density  consis¬ 
tent  with  the  solver  mesh. 

□  Field 

User  selects  a  calculation  type  from  a  pull-down  menu. 

■  SpeciesJVJntegral:  calculates  the  amount  of  species  that  exists  in  the  defined  volume  at  each  time. 
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Simulation  Results  -  Query 

When  the  query  is  complete,  a  Simulation  Results  window  opens  containing  both  NetFlow  and  Query  results. 


Closes  Visualizer  Control 
Panel  and  hierarchical 
results  window;  control  is 
returned  to  FlumeCAD 
Solvers  window. 


□  regionQuery  Table 

Opens  a  window  containing  a  table  of  results  from  the  In  a  Region  query.  Here,  regionQuery  is  the  label  (see 
Query  -In  a  Region  on  page  F2-142). 


□  regionQuery  Graph 

Opens  a  window  containing  graphical  results  from  the  In  a  Region  query.  Here,  regionQuery  is  the  label  (see 
Query  -  In  a  Region  on  page  F2-142). 
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Appendix  2:  Research  Publications 
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METROLOGY  AND  SIMULATION  OF  CHEMICAL  TRANSPORT  IN 

MICROCHANNELS 


P.  M.  St.  John*,  T.  Woudenberg,  and  C.  Connell 

PE  Applied  Biosystems 
850  Lincoln  Centre  Drive 
Foster  City,  CA  94404 


M.  Deshpande  and  J.  R.  Gilbert 

Microcosm  Technologies,  Inc. 

215  1st  St 

Cambridge,  MA  02142 

M.  Garguilo  and  P.  Paul 

Sandia  National  Laboratory 
Livermore,  CA 

ABSTRACT 

We  are  working  towards  building  a  CAD  tool  (NctFlow) 
for  microfluidic  systems.  To  support  the  development  of  this 
tool,  diffusion  and  flow  experiments  are  being  performed  in 
microchannels.  In  this  paper,  we  describe  the  novel  flow 
visualization  approach  developed  for  these  experiments  and 
report  experimental  results  and  corresponding  simulations  for 
fluid  transport  in  microchannels.  The  experimental  approach 
relies  on  the  use  of  a  caged,  charged  fluorescent  dye  which  can 
be  activated  (uncaged)  by  a  UV  laser  beam  and  tracked  by 
fluorescence  imaging  [1],  In  contrast  to  standard  methods  for 
inserting  a  dye  plug  in  a  microchannel,  this  is  the  only  technique 
which  allows  definition  of  very  precise  regions  at  any  place 
along  the  channel.  Specifically,  this  allows  a  narrow  definition 
of  starting  plugs  and  therefore  an  improved  determination  of  the 
flow  profile.  We  have  used  this  technique  to  study  diffusion, 
pressure-driven  flow,  and  electrokinetic  flow  in  various 
microchannel  geometries. 

INTRODUCTION 

There  is  a  wide  interest  in  micron-scale  integrated 
chemical/biochemical  analysis  or  synthesis  systems,  also 
referred  to  as  lab-on-a-chip.  The  microfluidic  systems  being 
explored  today  rely  on  both  pressure  [2]  and  electrokinetic 
effects  [3,4]  for  controlling  fluid  transport.  Complicated 
relationships  between  the  shape  and  size  of  the  microchannels 
and  the  behavior  of  multi-component  fluids  are  not  completely 
understood.  Researchers  are  forced  to  use  costly  trial  and  error 
methods  in  the  design  of  microfluidic  systems. 

Discrepancies  between  experimental  measurements  and 
simulations  exist  because  the  modeling  equations  may  not  fully 
capture  the  essential  physics  and  because  the  experimental 
parameters  used  in  the  simulations  are  often  approximate  (e.g. 
calculation  of  diffusion  constant  and  mobility  in  a 
microchannel).  We  are  attempting  to  simulate  simple  flow  in 
specific  microchannel  geometries  and  extract  parameters  which 
will  enable  us  to  predict  flow  in  more  complex  structures  and 
ultimately  allow  us  to  take  an  active  role  in  the  design  of 
microfluidic  structures.  Flow  predictions  have  direct 
applications  in  microcapillary  electrophoresis  systems  where 
long  channels  are  advantageous  for  resolved  separation  of  DNA 
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fragments  out  to  several  kilobases  but  device  area  must  be 
minimized.  Flow  predictions  arc  also  useful  for  determining  the 
geometries  needed  for  chemical  reactions  in  microchannels. 

EXPERIMENTAL 

Silicon  microchannel  structures  were  fabricated  as 
“negative”  masters  using  standard  photolithography  and  an  STS 
deep  reactive  ion  etch  (DRIE)  was  used  to  create  raised 
rectangular  structures  with  dimensions  of  50  pm  in  width,  5  cm 
in  length  and  40  -  80  pm  in  height.  An  RTV  silicone  elastomer 
(Dow  Coming)  was  used  to  form  positive  replicas  of  the  silicon 
structures  (Fig.  1)  [5].  Scanning  electron  micrographs  (SEMs) 
show  the  etched  silicon  sidewalls  and  the  corresponding 
elastomer  replica  (Au  sputtered)  (Fig.  2).  The  striae  located  near 
the  top  of  the  silicon  structure  was  a  result  of  over-etching.  The 
etched  silicon  region  was  smooth  to  an  RMS  roughness  of  ~ 
7.26  nm,  determined  by  atomic  force  microscopy.  The 
elastomer  SEM  shows  the  reproduction  of  the  features  present  in 


Figure  1.  Photograph  of  the  elastomer  device  with  three  test 
microchannels.  The  circles  are  the  injection  wells  and  the 
squares  are  used  to  study  pressure  driven  flow  in  a  sudden 
expansion. 
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Figure  2a.  SEM  of  the  silicon  mold 


Figure  2b:  SEM  of  the  elastomer  channel 


the  silicon.  Note  the  striations  which  run  the  length  of  the 
channel  on  the  lower  portion  of  the  channel  wall.  This 
corresponds  to  the  rough  overetched  region  of  the  silicon  mold. 
The  region  of  the  elastomer  which  replicated  the  etched  silicon 


surface  had  an  RMS  roughness  of  ~  7  nm.  The  rough  edge  at 
the  opening  of  the  elastomer  channel  came  from  slicing  the 
elastomer  structure  with  a  razor  blade.  The  channel  was  not 
cleaned  prior  to  metal  sputtering  and  particulates  can  be  seen  on 
the  channel  floor.  The  elastomer  channels  were  sealed  with 
either  glass  microscope  slides  or  elastomer  coated  glass 
microscope  slides  and  wells  were  created  in  the  elastomer  using 
a  mechanical  hole  punch  (Figure  1). 

Serpentine  microchannels  in  borofloat  glass  (not  shown) 
were  fabricated  by  Alberta  Microelectronics  Centre  (Alberta, 
Canada)  using  standard  photolithography  and  an  isotropic  wet 
etch  to  yield  a  smooth  hemispherical  cross-section.  A 
commercially  available  additive  made  from  a  polyacrylamide 
short  chain  polymer  (ABI  DNA  Fragment  Analysis  Reagent  7% 
w/w)  was  used  (0.1%  v/v  aqueous)  to  suppress  electroosmotic 
flow  (EOF)  [6]  in  both  glass  and  elastomer  structures.  Although 
EOF  suppression  was  verified  by  observing  ion  travel  under  an 
applied  field,  the  extent  of  suppression  was  not  quantitated.  In 
this  paper,  electrophoretic  flow  refers  to  negatively  charged  dye 
migration  towards  the  positive  electrode  and  electroosmotic 
flow  refers  to  dye  carried  toward  the  negative  electrode,  despite 
charge. 

The  caged  fluorophore  (fluorescein)  was  purchased  from 
Molecular  Probes  and  used  in  micromolar  quantities  dissolved  in 
either  distilled  water  or  Tris-EDTA  buffer,  pH  =  7.4.  The  dye 
is  negatively  charged,  both  caged  and  uncaged  (fluorescein,  -2). 
When  caged,  the  dye  solution  is  non-fluorescent.  Upon  photo¬ 
activation  of  a  UV  light  illuminated  volume,  the  protecting 
group  is  cleaved  and  the  dye  fluoresces.  Approximately  0.3  mW 
average  power  of  355  nm  light  from  a  pulsed  Nd:YAG  laser 
(Uniphase)  was  used  for  ~  0.1  s  to  uncage  the  fluorophore  in  a 
20  pm  spot.  In  the  electrokinetic  experiments,  the  laser  beam 
was  focused  using  a  UV  microscope  objective  which  defined  a 
sharp  start  zone.  A  continuous  wave  Microblue  diode  pumped 
laser  at  473  nm  (Uniphase)  was  used  for  fluorescence  excitation 
(Figure  3).  The  Microblue  was  positioned  close  to  Brewster’s 
angle  for  glass  to  minimize  laser  scatter  into  the  microscope. 
Pressure  driven  flow  was  imaged  at  Sandia  National 
Laboratories  where  the  experimental  setup  is  similar  to  that 
described  above  with  the  exceptions  that  the  uncaged  region  is 
rectangular  in  shape  and  that  the  uncaging  time  is  much  less 
than  one  millisecond.  Fluorescence  images  of  the  molecules 
diffusing  and  moving  in  the  local  flow  were  collected  using  a 
microscope  objective  (total  magnification  =  1 00  X)  and  a  video 
rate,  interlaced  camera.  Images  were  analyzed  and  flow 
parameters  were  extracted  for  simulations.  The  evolution  of  the 
fluorescent  profiles  under  electrokinetic  and  pressure-driven 
flow  are  shown  here  by  comparing  a  sequence  of  timed  images 
to  corresponding  times  from  the  simulations. 

The  numerical  simulations  were  performed  using  the 
NetFlow  module  in  Memcad  [7].  It  employs  a  three-dimensional 
finite  element  based  tool  to  solve  the  Navier-Stokes  equations. 
The  electrokinetic  effects  can  be  broadly  divided  into  two 
categories:  electrophoresis,  involving  the  transport  of  a  charged 
species  through  a  carrier  due  to  a  differential  in  the  mobility,  and 
electroosmosis,  involving  the  pumping  of  the  carrier  fluid  due  to 
wall  based  effects.  Electrokinetic  effects  were  incorporated  by 
coupling  the  Navier-Stokes  equations  with  the  Poisson  equation 
for  the  electric  field  and  with  the  Poisson-Boltzmann  equation 
for  the  zeta  potential  in  the  electroosmotic  cases  [8]. 

The  equations  are  solved  here  under  the  assumption  of  a 
dilute  solution  i.e.  the  carried  species  does  not  affect  the  material 
properties  of  the  carrier.  This  allows  the  equations  to  be  solved 
in  their  incompressible  form.  The  finite  element  meshes  were 
chosen  generally  to  allow  for  adequate  resolution  of  the 
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boundary  layer  and  accurate  representation  of  all  the 
experimentally  observed  fluid  physics. 

In  the  design  of  practical  micro-devices,  however, 
accurate  coupling  of  the  elcctrokinetic  effects  requires  the 
knowledge  of  several  material  and  physical  parameters,  such  as 
the  diffusivity  and  the  mobility  of  the  solute  in  the  solvent  and 
extent  of  surface  charge  on  the  inner  walls  of  the  microchanncls. 
Additionally  wall-based  effects  such  as  chemical  binding  sites 
might  exist  that  affect  the  flow  patterns.  Our  approach  is  to 
extract  the  required  material  parameters  through  quantitative 
comparison  with  experimental  measurements  in  simple 
geometries  and  use  the  extracted  parameters  for  simulations  and 
comparisons  in  more  complex  geometries. 


Schematic  of  the  uncaging  apparatus 


UV  pulse  initiates  dye 

I - ~1 

_ H _ L.,.*-- _ 

time  evolved 

Figure  3.  Top:  Schematic  of  the  experiment.  The  Nd.YAG 
laser  is  used  to  uncage  the  dye  and  a  Microblue  473  nm  laser  is 
used  to  excite  the  fluorescence.  Bottom:  Illustration  of  the 
uncaging  process.  The  starting  fluorescent  fluid  volume  is 
defined  by  the  uncaging  laser  (either  as  a  spot  or  as  a  ribbon). 

RESULTS  AND  DISCUSSION 

Figure  4  shows  a  timed  series  of  pressure-driven  flow  in 
microchannels  which  expand  into  a  large  planar  reservoir.  This 
figure  demonstrates  that  the  simulation  qualitatively  reproduces 
the  dye  dispersion  from  laminar  flow  seen  during  the 
experiment.  The  channel  width  was  100  pm  by  40  pm  deep. 
The  flow  rate  was  approximately  1 5  pL/hr. 

The  Reynolds  number.  Re,  in  this  experiment  is  typical  for 
bioanalytical  systems  (Re  =  0.05)  although,  it  appears  low  when 
compared  to  typical  numbers  in  micropumps  (Re  ~  100).  The 
images  therefore  differ  from  what  on  might  expect  to  see  at 
higher  Reynolds  number.  The  governing  equations  become 
linear  and  harmonic  which  results  in  the  flow  pattern  appearing 
as  if  it  evolved  from  a  point  source  in  Fig.  4.  The  experiments 
and  simulations  in  the  figure  show  good  agreement  which  is  to 
be  expected  for  low  Reynolds  number  flows.  Similar 
experiments  and  stimulations  were  conducted  on  various  other 
geometrical  microchannel  shapes  and  serve  to  verify  the 


capabilities  of  our  module  in  the  simulation  of  pressure-driven 

flows. 


Figure  4.  Experimental  and  simulated  pressure-driven  flow 
entering  a  sudden  expansion  in  an  elastomer  structure. 

Electrokinctic  measurements  using  a  range  of  applied 
fields  and  in  a  number  of  geometries  were  compared  to 
simulations  of  electroosmotic  and  electrophoretic  flow  in  both 
elastomer  and  glass  microchanncls.  Figure  5  shows  the  simple 
case  of  a  time-lapse  experimental  electrophoretic  flow  at  a  field 
strength  (E)  =  85.8  V/cm  in  a  50  fim  straight  elastomer  channel 
with  a  0.1%  flowable  and  water  soluble  polymer  added  to  the 
dye  solution  to  suppress  EOF  |6).  For  the  case  of 
electrophoretic  flow  where  EOF  is  assumed  negligible,  the 
simulated  and  experimental  band  shape  are  similar.  The 
simulated  image  (right)  shows  some  of  the  features  of  the 
NctFlow  software.  This  image  is  a  2-D  histogram  of  the 
integrated  experimental  fluorescent  profile  for  t  =  0.99  s.  The  2- 
D  histograms  provide  crucial  information  regarding  the  surface 
forces  acting  on  the  flow.  The  simulation  shows  an  intense 
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t  =  2.31 

Figure  5.  Left:  Experimental  electrophoretic  flow  in  an 
elastomer  channel,  E  -  85.8  V/cm.  Right:  Column  integration 
tool  from  NetFlow  used  to  form  2-D  image/histogram 
corresponding  to  the  fluorescence. 
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5.21  x  10"5  cm2/s  (at  298  K)  [9]  (sucrose  has  60%  the  mass  of 
fluorescein)  and  Dnuoresce;n  was  reported  to  be  3.3  x  10‘6  cm2/s  by 
polarography  (experimental  conditions  unknown)  |4].  The 
velocity  of  the  uncaged  dye  was  extracted  from  the  experimental 
data  by  plotting  the  distance  versus  time  for  a  fluorescent  plug 
moving  in  an  applied  field  (extracted  from  experimental  images 
like  those  shown  in  Fig.  5).  A  linear  regression  was  used  to 
obtain  the  slope  or  velocity  of  the  line.  Velocities  at  several 
voltages  were  obtained  in  this  way  and  plotted  versus  the  field 
strength.  In  the  absence  of  electroosmotic  flow,  the  slope  of  the 
line  is  equal  to  the  electrophoretic  mobility  (|i).  The  mobility 
was  estimated  to  be  2.3  x  10'4cm2/Vs,  using  the  linear  region  of 
the  graph  (E  >  60  V/cm).  The  plot  shows  a  non-linear  behavior 
between  velocity  and  field  strength  for  E  <  60  V/cm.  We 
suspect  this  is  not  an  interaction  between  the  dye  molecules  and 
the  solvent  because  the  dye  is  in  such  low  concentration 
compared  to  the  solvent.  A  more  likely  explanation  might  be 
the  error  associated  with  obtaining  measurements  at  low  field 
strengths  when  the  diffusion  velocity  is  approximately  the  same 
order  of  magnitude  as  the  convective  velocity.  The  extracted 
diffusivity  and  mobility  from  the  elastomer  straight  channels 
was  used  successfully  in  modeling  the  plug  profile  in  the  glass 
u-shaped  channels  implying  that  these  parameters  can  be  used  to 
predict  flow  in  a  variety  of  structures. 
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Figure  7.  A  plot  of  the  velocity  vs.  field  strength  used  to  obtain 
the  mobility  of  uncaged  fluorescein. 


Diffusivity  (cmA2/s) 

Figure  6.  Graph  of  the  simulated  band  width  at  two  distinct 
times  (squares  and  triangles)  vs.  a  range  of  dijfusivities. 
Experimental  data  (circles)  were  added  to  the  plot  to  extract  the 
diffusivity 

center  region  (high  mass  fraction)  and  regions  of  very  low 
intensity  along  the  channel  walls,  implying  that  the  fluorescent 
plug  may  be  repelled  in  this  region.  The  band  broadening 
occurs  rapidly  because  the  low  mass  dye  has  a  high  diffusion 
constant,  ~  1.3  x  10'6  cm2/s.  The  migration  of  the  band  is 
towards  the  positive  electrode  which  is  expected  for 
electrophoretic  flow  of  a  negatively  charged  species  with 
suppressed  EOF. 

The  simulations  can  also  be  used  to  extract  parameters 
from  the  experimental  data.  Figure  6  shows  a  graph  of  the 
simulated  band  width  at  two  distinct  times  (t  =  0.99  s,  triangles, 
and  t  =  1 .98  s,  squares)  as  a  function  of  predicted  diffusivities. 
We  were  able  to  extract  a  diffusivity  of  the  uncaged  fluorophore 
in  the  microchannel,  D  =  1.3  x  10"6  cm2/s,  by  comparing  the 
experimental  bandwidths  (circles)  at  t  =  1.98  s  and  t  =  0.99  s 
with  the  calculated  curves.  This  number  seems  reasonable 
considering  other  reported  numbers  for  diffusivity,  Dsucrose  = 


Although  the  field  strengths  used  here  are  typical  for  DNA 
electrophoresis  (-100  V/cm),  the  diffusivity  of  the  caged 
fluorophore  in  aqueous  solution  is  much  higher  than  that  of 
DNA  because  of  its  lower  molecular  weight  (m.w.caged  fluorescein  = 
827  d,  m.w.  DNA  -  several  hundred  kd).  Therefore,  the  band 
profile  is  dominated  by  diffusion  at  these  field  strengths.  In 
order  to  predict  the  loss  of  resolution  from  diffusion  in  DNA 
separation,  a  lower  diffusivity  must  be  used.  We  have  also 
examined  the  flow  profile  in  more  complex  geometries.  Figure 
8  compares  the  experimental  to  simulated  fluorescence  profile 
under  electrophoresis  in  a  u-shaped  glass  microchannel  at  four 
different  time  steps.  The  images  corresponds  to  E  =  69.5  V/cm. 
The  inset  shows  the  current  density  in  the  same  section  of  the 
microchannel.  As  the  figure  shows,  the  higher  current  density 
on  the  inner  wall  of  the  bend  causes  the  inner  region  of  the  band 
to  move  faster  than  the  outer  region,  resulting  in  an  asymmetric 
profile  and  eventual  band  broadening.  The  greater  distance 
traveled  by  the  species  along  the  outer  wall  further  adds 
broadening  and  asymmetry  to  the  band.  As  it  travels  past  the  u- 
shaped  bend  the  band  takes  the  shape  of  a  parallelogram.  In 
principle,  this  shape  can  be  corrected  by  allowing  the  band  to 
negotiate  an  equal  and  opposite  turn  downstream  of  this  bend. 
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Figure  8.  Experimental  (left)  and  simulated  (right) 
electrophoretic  flow  in  a  50  pm  diameter  glass  microchannel 
(radius  of  curvature  ~  0.25  mm),  E  =  69.5  V/cm  (inset:  current 
density  in  channel,  the  lightest  region  corresponds  to  the  highest 
current  density). 

However,  the  broadening  effects  of  the  bend  and  the  diffusion 
cannot  be  corrected. 

The  simulation  at  E  =  69.5  V/cm  in  Fig.  8  closely  matched 
the  experiment.  As  the  field  strength  was  increased  further  the 
band  was  observed  to  take  on  a  parabolic  shape  that  grew 
increasingly  well-defined  with  increasing  field  strength.  This 
behavior  is  contrary  to  expectations  and  is  not  represented  in  the 
numerical  models.  Consequently,  the  agreement  between 
experiment  and  simulations  falls  off  at  high  voltages.  Figure  9 
shows  the  experimental  band  profile  at  0.66  s  after  uncaging  the 
fluorophore  at  a  range  of  field  strengths,  E  =  1 08  -  1 25  V/cm. 

As  the  field  strength  increases  above  E  ~  100  V/cm,  the 
parabolic  profile  becomes  more  pronounced.  The  increasing 
parabolic  profile  with  increasing  field  strength  has  also  been 
observed  in  straight  microchannels  (data  not  shown)  and  straight 
capillaries  [10]  and  is  therefore  not  due  to  the  asymmetry  in  the 
field.  The  simulations  in  this  field  range  did  not  match  the 


experimental  data.  There  arc  several  possible  reasons  for  the 
observed  parabolic  shape  at  high  voltages.  One  possibility  is  the 
variation  in  properties  (diffusivity,  viscosity,  and  mobility) 
across  the  channel  due  to  joule  heating  which  would  cause  a 
temperature  gradient  in  the  channel,  that  becomes  more 
significant  at  higher  voltages.  Temperature  gradients  and  their 
effects  on  the  electrokinetic  flow  behavior  in  microchanncls 
have  been  observed  by  other  researchers  [11).  A  1  -  2°  C  rise  in 
temperature  in  the  channel  can  cause  the  mobility  (and  viscosity 


Figure  9  A  scries  of  profiles  at  l  =  0.66  s  after  the  dye  was 
uncaged  at  varying  field  strengths. 

of  the  carrier)  to  change  by  2  -  3%.  At  higher  voltages,  where 
the  temperature  rise  is  significant,  the  change  in  these  material 
properties  may  affect  the  overall  flow  pattern.  A  temperature 
rise  may  also  effect  the  polymer  additive  [6],  resulting  in  an 
incomplete  EOF  suppression.  If  EOF  is  present,  then  non- 
uniformities  in  the  zeta  potential  could  induce  a  pressure 
gradient  and  a  corresponding  parabolic  component  of  the  flow. 
Finally,  the  parabolic  shape  may  be  present  at  all  field  strengths 
but  masked  by  the  more  rapid  decrease  in  fluorescence  at  lower 
field  strengths  due  to  diffusion.  We  arc  currently  exploring 
these  mechanisms  to  explain  the  observed  deviations  from  plug¬ 
like  flow,  using  both  experimental  and  numerical  techniques. 
An  appropriate  physical  model  can  then  be  developed  and 
incorporated  into  the  electrokinetic  terms  to  represent  this 
behavior. 

CONCLUSION 

Wc  have  used  a  caged  fluorophore  to  image  electrokinetic 
and  pressure-driven  flow  profiles  in  both  glass  and  elastomer 
microchannels  with  varying  geometries.  Parameters  such  as 
diffusivity  and  mobility  were  extracted  from  the  experimental 
data  and  used  in  the  simulations.  This  imaging  technique  has 
already  revealed  some  interesting  and  unexpected  behavior  in 
these  relatively  simple  fluidic  systems.  We  are  looking  forward 
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to  using  this  method  to  study  these  behaviors  and  extract 
theoretical  models  or  behavioral  models  suitable  for  insertion  in 
modeling  tools  such  as  NetFlow. 
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ABSTRACT 

This  paper  presents  a  numerical  framework  for  design-based  analyses  of  electrokinetic  flow  in  interconnects. 
Electrokinetic  effects,  which  can  be  broadly  divided  into  electrophoresis  and  electroosmosis,  are  of  importance  in  providing 
a  transport  mechanism  in  microfluidic  devices  for  both  pumping  and  separation.  Models  for  the  electrokinetic  effects  can  be 
derived  and  coupled  to  the  fluid  dynamic  equations  through  appropriate  source  terms.  In  the  design  of  practical  micro¬ 
devices,  however,  accurate  coupling  of  the  electrokinetic  effects  requires  the  knowledge  of  several  material  and  physical 
parameters,  such  as  the  diffusivity  and  the  mobility  of  the  solute  in  the  solvent.  Additionally  wall-based  effects  such  as 
chemical  binding  sites  might  exist  that  affect  the  flow  patterns. 

In  this  paper,  we  address  some  of  these  issues  by  describing  a  synergistic  numerical/experimental  process  to  extract  the 
parameters  required.  Experiments  were  conducted  to  provide  the  numerical  simulations  with  a  mechanism  to  extract  these 
parameters  based  on  quantitative  comparisons  with  each  other.  These  parameters  were  then  applied  in  predicting  further 
experiments  to  validate  the  process  As  part  of  this  research,  we  have  created  NetFlow,  a  tool  for  micro-fluid  analyses.  The 
tool  can  be  validated  and  applied  in  existing  technologies  by  first  creating  test  structures  to  extract  representations  of  the 
physical  phenomena  in  the  device,  and  then  applying  them  in  the  design  analyses  to  predict  correct  behavior. 

Keywords  :  Electrophoresis,  Electroosmosis  ,  Numerical  Models,  Design  Analysis,  CAD. 


1.  INTRODUCTION 

There  is  a  wide  interest  in  micron-scale  integrated  chemical/biochemical  analysis  or  synthesis  systems,  also  referred  to 
as  lab-on-a-chip.  The  microfluidic  systems  being  explored  today  rely  on  both  pressure1  and  electrokinetic  effects2'4  for 
controlling  fluid  transport.  The  design  of  such  systems  is  hampered  by  the  lack  of  information  required  to  accurately 
simulate  the  systems.  The  problem  is  further  complicated  by  the  differing  solvent/solute  combinations  that  are  employed  in 
different  applications  as  well  as  the  actual  material  and  manufacturing  process  of  the  system  components.  Design  based 
analyses  of  these  systems  therefore  require  a  priori  extraction  of  the  material  properties.  In  many  cases,  however,  the 
material  properties,  such  as  the  diffusion  constant  are  difficult  to  extract  using  purely  experimental  techniques.  On  the  other 
hand,  a  combined  experimental/simulation  technique  is  ideally  suited  for  this  purpose,  and  is  the  focus  of  this  paper. 


1  Further  author  Information  - 

Manish  Deshpande  (correspondence):  manish@memcad.com.  215  First  Street,  Cambridge,  MA  02142. 
Tel  (617)  225-0094  x222,  www.memcad.com  for  information. 
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Simulation  of  electrokinetic  flows  have  been  reported  by  several  researchers  in  the  literature  for  both  electrophoretic5,6 
and  electroosmotic'  flows.  In  most  cases,  these  simulations  have  focussed  on  the  numerical  problem  with  well-defined 
material  properties.  Using  the  treatment  of  electrokinetic  flow  presented  there  as  a  baseline,  we  demonstrate  our  attempt  at 
the  inverse  problem  -  i.e.,  the  extraction  of  the  parameters  required  to  adequately  replicate  the  experiment.  With  these 
parameters  in  hand,  we  can  then  proceed  to  more  complicated  design  questions  that  can  be  answered  by  simulations.  In 
effect,  we  demonstrate  a  mechanism  to  extract  the  reduced-order  component  model  for  microchannel  flow  that  can  then  be 
inserted  into  the  full  system  model. 


2  GOVERNING  EQUATIONS 

The  basic  equations  describing  the  fluid  motion  are  the  Navier-Stokes  equations  with  appropriate  electromigratory  flux 
terms  to  represent  the  effect  of  the  applied  electric  field  on  the  carrier  and/or  the  charged  species.  The  effect  of  the  applied 
field  can  be  divided  into  two  fundamental  components  : 

Electrophoresis  :  The  basis  for  electrophoresis  is  the  differential  migration  of  the  charged  species  ions  relative  to  the 
carrier  molecules  under  the  application  of  the  external  field.  The  differential  migration  is  primarily  an  effect  of  the 
difference  in  the  net  charge  between  the  solvent  and  solute  ions,  although  frictional  effects  may  also  have  some  relevance. 
The  migration  velocity  of  the  charged  species  can  be  expressed  in  terms  of  the  applied  field  strength  as 

%  =  **  £  E>] 

where  pcp  is  the  electrophoretic  mobility  of  the  ion  in  the  carrier  species.  It  is  important  to  note  that  in  most  cases  the 
carrier  does  not  move  under  electrophoresis. 

Electroosmosis  :  Electroosmosis,  in  contrast,  is  a  macroscopic  phenomenon  involving  the  pumping  of  a  fluid  through  a 
channel  under  the  application  of  the  field.  In  most  cases,  walls  in  microchannels  are  characterized  by  the  presence  of 
surface  charges.  The  charge  may  either  be  due  to  the  property  of  the  wall  or  by  adsorption  of  the  charged  species  from  the 
buffer.  In  the  presence  of  an  electrolyte  the  surface  charge  density  induces  the  formation  of  a  double  layer  in  the  fluid  by 
attracting  oppositely  charged  ions  from  the  electrolyte  to  the  immediate  vicinity  of  the  wall.  The  application  of  the  electric 
field  exerts  a  force  on  the  fluid,  that  is  initially  felt  only  within  the  double  layer.  As  a  result  the  fluid  in  the  near  vicinity  of 
the  wall  starts  to  move.  Due  to  the  viscous  forces  the  fluid  in  the  center  of  the  channel  is  also  accelerated  until  the  net 
velocity  gradient  in  the  radial  direction  is  zero  and  the  whole  fluid  in  the  channel  moves  at  a  constant  velocity. 

The  determination  of  the  electroosmotic  flow  field  requires  the  solution  of  the  Navier-Stokes  equation.  Incorporating 
the  electroosmotic  effect  as  a  force  the  equations  become 

V  •  V  =0  (Continuity) 
eo 

D(pV  ) 

-  =  -  Vp  +  V»  /rV  V  +  p  E  [2] 

Dt  eo  e 

Where  Veo  is  the  induced  electroosmotic  velocity.  D/Dl  is  the  material  derivative  of  the  momentum,  (a  is  the  fluid 
viscosity  and  p  is  the  pressure. 

The  last  term  in  the  momentum  equation  represents  the  electroosmotic  force  on  the  fluid.  Here  pc  is  the  charge  density 
and  E  is  the  electric  field  intensity.  The  electric  field  can  be  determined  by  the  solution  of  the  potential  equation 

£  =  -V$  V2<h  =  ~(pe/s  )  [3] 

The  electric  potential  <f>  can  be  further  decomposed  into  two  components,  (p.  the  external  applied  potential  and  the 
zeta  potential  at  the  walls. 


=  <j>  +  £ 


[4] 
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Under  the  assumption  that  the  zeta  potential  effects  are  confined  to  a  very  small  region  near  the  wall,  the  charge 
distribution  can  be  assumed  to  be  governed  by  the  zeta  potential  alone,  independent  of  the  external  field.  This  allows  the 
decomposition  of  Eq.  [3]  into  separate  equations  for  the  applied  field  and  the  zeta  potential 

V2^  =0  (Applied  Potential)  [5a] 

VV  =  -( pe!s )  (Zeta  Potential)  [5£>] 


The  equation  for  the  zeta  potential  can  be  simplified  by  the  Debye-Hiickel8  treatment  for  the  charge  density.  This  is 
based  on  the  concept  of  a  diffuse  double  layer,  proposed  initially  by  Gouy  and  Chapman8,  which  assumes  that  the  double 
layer  extends  for  some  finite  distance  into  the  fluid.  Through  the  Debye-Hiickel  approximation  the  charge  density  on  the 
walls  can  be  determined  through  the  Boltzmann  equation  and  results  in  the  Poisson-Boltzmann  equation  for  the  zeta 
potential 
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where  n0  is  the  ionic  concentration,  z  is  the  valence  of  the  charged  buffer,  e  is  the  electron  charge,  k  is  the  Boltzmann 
constant  and  T  is  the  temperature.  Eo  is  the  permittivity  of  vacuum. 


Under  further  simplification  under  the  assumption  that  ze  «  kT  equation  [6]  reduces  to  the  form 


£  =  Kz  £ 


[7] 


where 


K  =  (2n0zV  lss0kTf 12 

is  known  as  the  Debye-Hiickel  parameter  and  “1/K”  is  termed  as  the  characteristic  thickness  of  the  double  layer. 

The  coupled  solution  of  the  Eqs  (3),  (5a)  and  (7)  yields  the  flow  field  under  the  conditions  specified. 

The  motion  of  a  charged  species  in  the  electric  field  can  can  be  determined  by  coupling  the  species  mass  equation  (8) 
with  the  above  equations. 


—  +fv  +V  \Cj 

dt  {  eo  eP  J  1 


DV2  Cj 
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The  convective  transport  of  the  species  is  through  the  combined  effect  of  the  electroosmotic  motion  of  the  carrier  fluid 
and  the  electrophoretic  transport  of  the  species  under  the  effect  of  the  applied  electric  field. 

The  equations  presented  above  are  derived  under  the  following  assumptions  - 

Dilute  Carrier  :  The  carrier  fluid  is  assumed  to  be  electroneutral  everywhere,  except  within  the  double  layer.  The 
charge  distribution  is  therefore  confined  to  a  small  region  within  the  double  layer. 

Dilute  Solution  :  The  carrier  fluid  is  assumed  to  be  the  predominant  species  in  the  mixture.  That  is,  the  mass  fraction  of 
the  charged  species  is  assumed  negligible  in  comparison  to  the  fraction  of  carrier.  This  is  generally  valid  in  most 
electrokinetic  problems  where  the  species  under  observation  are  generally  in  the  millimolar  range  (corresponding  to  mass 
fractions  of  10'3  and  lower).  This  implication  of  the  dilute  assumption  is  that  the  solute  species  does  not  affect  the  material 
properties  of  the  mixture.  This  allows  the  species  transport  equation  to  be  decoupled  from  the  momentum  equation. 
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Individual  Species  do  not  affect  each  other  :  This  assumption  is  used  in  describing  the  flux  terms  in  the  above  equations, 
in  that  the  flux  of  one  species  does  not  depend  on  the  fluxes  of  the  other  species.  This  is  only  true  in  dilute  solutions. 

No  Chemical  Reactions  :  The  charged  species  are  assumed  to  be  fully  ionized  in  the  mixture. 

As  mentioned  above  the  dilute  solution  allows  the  density  of  the  mixture  to  be  assumed  constant  and  equal  to  the 
density  of  the  carrier.  This  reduces  the  problem  to  the  incompressible  form.  The  decoupling  of  the  species  equation  from 
the  momentum  equation  also  implies  that  the  species  conservation  condition  is  not  relevant  and  can  be  discarded. 


3.  THE  DESIGN  PROBLEM 


The  design  analysis  for  a  device  of  interest  generally  involves  the  representation  of  the  device  as  a  system  comprising 
several  individual  components.  The  physics  of  each  component  is  then  extracted  from  detailed  computation  for  that 
component.  The  physics  of  the  component  can  be  represented  in  the  system,  yielding  a  reduced  order  model  for  the  device. 

In  a  microfluidic  system  the  components  are  traditionally  active  and  passive.  Reaction  sites,  switching  chambers  etc. 
form  active  components,  whereas  transport  and  separation  channels  form  passive  components.  Netflow  is  about  the  design 
of  these  passive  components. 

The  numerical  framework  described  above  enables  the  computation  of  the  electric  and  fluid  flow  fields  for  a  passive 
component  under  a  specified  set  of  conditions.  A  design  analysis  for  a  practical  device,  on  the  other  hand,  presents  several 
unique  challenges.  These  can  be  both  numerical  and  physical.  Accurate  modeling  of  the  device  physics  requires  accurate 
material  property'  definitions.  In  general,  the  carrier  fluid  is  a  standard  buffer  with  several  additives  that  are  specific  to  the 
application  under  design.  The  material  properties  of  the  carrier  fluid  can  be  measured  relatively  easily.  The  properties  of  the 
solute  species  such  as  its  diffusivity  in  the  carrier  fluid  and  its  electrophoretic  mobility  are,  however,  considerably  harder  to 
extract.  In  addition,  these  properties  are  often  characteristic  to  specific  design  conditions. 

A  second  issue  is  the  effect  of  the  chemistry  of  the  channel  walls.  In  microchannels,  the  walls  of  the  channels  often  have 
non-uniform  chemistry  on  them  as  a  result  of  the  manufacturing  process,  through  additives  or  other  reasons.  This  non- 
uniform  wall  chemistry  results  in  flow  behavior  that  cannot  be  obtained  from  the  above  equations  without  the  addition  of  a 
model  to  represent  the  effect.  For  example,  the  presence  of  chemical  binding  sites  on  walls  can  result  in  the  retardation  of 
the  species  in  the  vicinity  of  the  wall.  This  retardation  is  highly  undesirable  since  it  leads  to  band  broadening, 
contamination  and  eventual  degradation  in  the  resolution  of  the  detection  signal.  The  effects  of  the  wall  chemistry  can  be 
represented  in  a  numerical  model  through  retardation  terms  in  the  appropriate  equations.  For  example,  the  affinity  between 
the  chemical  present  on  the  wall  and  the  charged  species  can  be  represented  as  a  retardation  force  on  the  charged  species. 
The  convection  of  the  charged  species  is  therefore  retarded  and  can  be  modeled  as  a  variation  in  the  convective  velocity 
that  is  proportional  to  its  distance  from  the  wall.  The  exact  dependence,  however,  can  only  be  extracted  by  quantitative 
comparison  between  experiment  and  simulation.  The  numerical  framework  presented  in  this  paper  enables  extraction  of  this 
dependency  as  well. 

In  the  case  of  electroosmosis,  the  accurate  solution  of  the  Poisson-Boltzmann  equation  requires  the  resolution  of  the 
double  layer.  The  thickness  of  this  double  layer  is  as  described  in  eq.  [7]  above.  In  general,  the  double  layer  thickness  is 
extremely  small  in  comparison  to  the  channel  width.  For  example,  under  nominal  conditions  for  a  50  micron  channel 
containing  water,  the  double  layer  thickness  is  of  O(nm).  Adequate  resolution  of  scales  of  this  order  are  not  generally 
tractable  and  further  simplification  is  necessary.  The  effect  of  the  potential  across  the  double  layer  (the  zeta  potential)  is  felt 
as  a  force  on  the  fluid  at  the  edge  of  the  double  layer.  The  solution  to  the  Poisson-Boltzmann  equation  can  then  be  reverted 
to  the  Smoluchowski  equation  at  the  outer  edge  of  the  double  layer : 


V 

eo 


—  Aeo  E 


[9] 


Which  is  analogous  to  the  equation  for  the  electrophoretic  velocity,  Eq.  [1],  The  electroosmotic  mobility  can  then  also 
be  extracted  using  an  approach  similar  to  the  electrophoretic  case. 

In  summary,  the  design  problem  is  as  follows  :  For  a  given  application,  the  various  physical  properties  that  affect  the 
behavior  of  the  device  must  be  extracted  in  an  appropriate  parameter  space.  The  properties  can  then  be  used  in  that 
parameter  space  to  conduct  further  analyses.  Netflow  is  the  first  step  in  enabling  this  process. 
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4.  SOLUTION  METHODOLOGY 


The  modeling  of  electrokinetic  effects  is  incorporated  into  the  NetFlow  module  in  the  MEMCAD  system.  MEMCAD  is 
an  integrated  design  environment  consisting  of  3D  design,  modeling  and  simulation  software  tools  which  enable  the 
creation  of  complex  MEMS/MST  devices.  Inherent  in  the  design  process  in  Memcad  is  the  ability  to  translate  from  a  layout 
and  process  view  of  the  device  to  a  solid  model  and  a  suite  of  simulators  to  characterize  the  various  physical  phenomena 
present  in  the  device.  Netflow,  which  is  a  part  of  the  suite,  is  specific  to  the  analysis  of  chemical  transport  in  microchannels. 
It  is  an  ongoing  program,  and  is  a  result  of  a  collaboration  between  Microcosm,  Perkin-Elmer  Corp.  and  Stanford 
University. 

The  numerical  solution  of  the  above  equations  is  performed  using  a  three-dimensional  finite  element  based  CFD  engine. 
The  Memcad  interface  allows  the  user  to  define  the  geometry  from  layout  and  process,  or,  optionally  directly  build  the 
model  in  the  solid  modeler.  Once  meshed,  the  model  is  available  to  a  suite  of  solvers  for  various  physical  domains,  and  to 
the  simulation  manager  to  parametrize  the  dependency  in  any  domain.  A  schematic  of  the  solution  methodology  is  shown 
in  Fig.  1. 

The  specific  problems  presented  here  are  for  electrophoretic  flow  in  microchannels.  The  parameters  that  need  to  be 
extracted  are  the  diffiisivity  and  mobility  of  the  charged  species  in  the  channels.  The  assumption  here  is  that  the  mobility 
and  diffusivity  are  independent  of  the  geometry  of  the  channels,  and  in  the  electrophoretic  case,  also  independent  of  the 
channel  wall  material.  That  is,  properties  extracted  from  simple  structures  can  be  applied  to  problems  in  more  complex 
structures,  assuming  the  channel  materials  and  manufacturing  processes  remain  the  same. 

5.  EXPERIMENTAL  SETUP 

Experimental  studies  of  electrophoretic  flow  in  microchannels  were  conducted  at  the  Applied  Biosystems  Division  at 
Perkin-Elmer,  and  at  Stanford.  A  brief  description  is  provided  here.  Greater  details  can  be  found  elsewhere9.  The  elastomer 
microchannels  were  fabricated  by  making  positive  replicas  of  raised  silicon  structures  fabricated  using  photolithography 
and  an  STS  deep  reactive  ion  etch  (DRIE).  The  elastomer  channels  were  sealed  with  either  glass  microscope  slides  or 
elastomer  coated  glass  microscope  slides. 

A  commercially  available  additive  made  from  a  polyacrylamide  short  chain  polymer  (ABI  DNA  Fragment  Analysis 
Reagent  7%  w/w)  was  used  (0.1%  v/v  aqueous)  to  suppress  electroosmotic  flow10  (EOF)  in  both  glass  and  elastomer 
structures.  Although  EOF  suppression  was  verified  by  observing  ion  travel  under  an  applied  field,  the  extent  of  suppression 
was  not  quantified.  In  this  paper,  electrophoretic  flow  refers  to  negatively  charged  dye  migration  towards  the  positive 
electrode  and  electroosmotic  flow  refers  to  dye  carried  toward  the  negative  electrode,  despite  charge. 

The  species  tracked  was  a  caged  fluorophore  (fluorescein)  was  purchased  from  Molecular  Probes  and  used  in 
micromolar  quantities  dissolved  in  either  distilled  water  or  Tris-EDTA  buffer,  pH  =  7.4.  The  dye  is  negatively  charged, 
both  caged  and  uncaged  (fluorescein,  -2).  When  caged,  the  dye  solution  is  non-fluorescent.  Upon  photo-activation  of  a  UV 
light  illuminated  volume,  the  protecting  group  is  cleaved  and  the  dye  fluoresces.  Fluorescence  images  of  the  molecules 
diffusing  and  moving  in  the  local  flow  were  collected  using  a  microscope  objective  (total  magnification  =  100  X)  and  a 
video  rate,  interlaced  camera.  Images  were  analyzed  to  extract  flow  parameters  by  comparing  against  simulations. 


6.  RESULTS 

Electrokinetic  measurements  using  a  range  of  applied  fields  and  in  a  number  of  geometries  were  compared  to 
simulations  of  electroosmotic  and  electrophoretic  flow  in  both  elastomer  and  glass  microchannels.  Representative  results 
are  shown  in  this  section.  The  field  strengths  studied  ranged  from  0  to  150  V/cm.  Simulations  were  conducted  for  the  same 
conditions  and  the  results  were  compared  to  experiment  to  extract  the  required  material  parameters  -  in  this  case  the 
diffusivity  and  mobility  of  the  charged  species. 

The  zero  voltage  case  enables  the  extraction  of  diffusivity  data,  since  diffusion  is  the  only  transport  mechanism  present. 
Note  that,  as  mentioned  before,  a  polymer  is  added  to  suppress  electroosmotic  flow  -  the  bulk  carrier  is  therefore  stationary. 
The  mechanism  for  extraction  of  the  diffusion  constant  is  as  follows  -  a  series  of  runs  are  conducted  for  a  range  of  diffusion 
constants  using  the  Simulation  Manager.  The  results  are  compared  with  the  experimental  observations  -  in  this  case  the 
broadening  of  the  band  with  time  -  and  the  diffusion  constant  is  then  extracted  corresponding  to  the  value  that  matches 
experiment. 
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Time  sequence  experimental  images  of  flow  at  zero  voltage  are  shown  in  Fig.  2.  The  corresponding  simulations  for  one 
diffusion  constant  are  also  shown  in  the  figure.  The  initial  band  shape  in  the  experiment  is  approximately  Gaussian  -  the 
initial  species  distribution  in  the  simulations  were  also  specified  to  be  Gaussian  to  represent  the  experiment.  The  band 
broadens  rapidly  because  the  dye  has  a  high  diffusion  constant  and  the  initial  Gaussian  distribution  quickly  becomes 
uniform  and  rectangular.  As  expected,  the  simulations  also  replicate  this  phenomenon. 

Figure  3  shows  a  graph  of  the  simulated  band  width  at  two  distinct  times  (t  =  0.99  s  and  t  =  1 .98  s)  as  a  function  of 
predicted  diffusivities.  The  diffusion  constant  can  be  extracted  from  this  graph  by  superimposing  the  experimentally 
measured  band  widths  at  the  same  times  (shown  as  the  circles  in  the  figure).  The  diffusion  constant  extracted  is 
approximately  1.3  x  10'6  cm2/s.  As  the  figure  shows  the  diffusion  constant  extracted  at  one  time  is  corroborated  by  that 
extracted  at  a  different  time.  The  constant  extracted  seems  reasonable  considering  other  reported  numbers  for  diffusivity, 
Dsucrose  =  =  5.21  x  10‘5  cm2/s  (at  298  K)11  (sucrose  has  60%  the  mass  of  fluorescein)  and  Dn110resccin  was  reported  to  be  3.3  x 
10'6  cm2/s  by  polarography  (experimental  conditions  unknown)3.  A  point  to  note  is  that  the  curves  corresponding  to  the 
simulations  are  extracted  from  a  single  Memcad  run  without  user  intervention. 

To  extract  the  mobility,  we  chose  a  field  strength  (E)  of  85.8  V/cm  for  the  comparison.  Figure  4  shows  the  experimental 
time  lapse  images  for  this  field  strength  in  similar  channels  as  above.  For  the  case  of  electrophoretic  flow  where  EOF  is 
assumed  negligible,  the  simulated  and  experimental  band  shape  are  similar.  The  simulated  image  (right)  shows  some  of  the 
features  of  the  NetFlow  software.  This  image  is  a  2-D  histogram  of  the  integrated  experimental  fluorescent  profile  for  t  = 
0.99  s.  The  2-D  histograms  provide  crucial  information  regarding  the  surface  forces  acting  on  the  flow.  The  simulation 
shows  an  intense  center  region  (high  mass  fraction)  and  regions  of  very  low  intensity  along  the  channel  walls,  implying  that 
the  fluorescent  plug  may  be  repelled  in  this  region.  The  band  broadening  occurs  rapidly  because  the  low  mass  dye  has  a 
high  diffusion  constant,  ~  1.3  x  10 6  cmz/s.  The  migration  of  the  band  is  towards  the  positive  electrode  which  is  expected 
for  electrophoretic  flow  of  a  negatively  charged  species  with  suppressed  EOF. 

The  velocity  of  the  uncaged  dye  was  extracted  from  the  experimental  data  by  plotting  the  distance  versus  time  for  a 
fluorescent  plug  moving  in  an  applied  field  shown  in  Figure  5(extracted  from  experimental  images  like  those  shown  in  Fig. 
4).  A  linear  regression  was  used  to  obtain  the  slope  or  velocity  of  the  line.  Velocities  at  several  voltages  were  obtained  in 
this  way  and  plotted  versus  the  field  strength.  In  the  absence  of  electroosmotic  flow,  the  slope  of  the  line  is  equal  to  the 
electrophoretic  mobility  (u).  The  mobility  was  estimated  to  be  2.3  x  KT4  cm2/Vs,  using  the  linear  region  of  the  graph  (E  > 
60  V/cm).  The  plot  shows  a  non-linear  behavior  between  velocity  and  field  strength  for  E  <  60  V/cm.  We  suspect  this  is 
not  an  interaction  between  the  dye  molecules  and  the  solvent  because  the  dye  is  in  such  low  concentration  compared  to  the 
solvent.  A  more  likely  explanation  might  be  the  error  associated  with  obtaining  measurements  at  low  field  strengths  when 
the  diffusion  velocity  is  approximately  the  same  order  of  magnitude  as  the  convective  velocity. 

Although  the  field  strengths  used  here  are  typical  for  DNA  electrophoresis  (-100  V/cm),  the  diffusivity  of  the  caged 
fluorophore  in  aqueous  solution  is  much  higher  than  that  of  DNA  because  of  its  lower  molecular  weight  (m.w.c.,,,c(i  fluort.sccin 
=  827  d,  m.w.  DNA  ~  several  hundred  kd).  Therefore,  the  band  profile  is  dominated  by  diffusion  at  these  field  strengths.  In 
order  to  predict  the  loss  of  resolution  from  diffusion  in  DNA  separation,  a  lower  diffusivity  must  be  used. 

The  success  of  this  approach  is  however  dependent  on  the  ability  to  use  the  extracted  parameters  to  flow  in  more 
complex  geometries.  To  demonstrate  this  we  chose  u-shaped  microchannels  using  the  same  buffer  and  dye  combination. 
Figure  6  4  shows  a  schematic  (top)  of  the  various  geometries  used  in  these  flow  studies.  The  motivation  in  using  a 
serpentine  structure  for  DNA  separation  is  that  we  can  separate  longer  fragments  of  DNA  in  nearly  the  same  space  required 
to  separate  the  shorter  fragments.  The  channels  were  fabricated  by  etching  a  series  of  50  pm  width,  various  length 
serpentine  channels  into  a  4”  x  4”  borosilicate  glass  plate  which  was  bonded  to  a  second  plate  to  seal.  These  channels  are 
trapezoidal  in  cross-section  and  require  considerable  resolution  to  accurately  simulate  the  flow  field.  Full  3D  computations 
were  conducted.  A  nominal  mesh  used  in  the  simulations  contained  roughly  2000  parabolic  finite  elements.  Initial 
simulations  to  verify  mesh  independence  were  conducted  under  nominal  conditions  and  are  not  presented  here. 

Figure  7  compares  the  experimental  to  simulated  fluorescence  profile  under  electrophoresis  in  a  u-shaped  glass 
microchannel  at  four  different  time  steps.  The  images  corresponds  to  E  =  69.5  V/cm.  As  before,  EOF  was  suppressed  by 
the  polymer  additive.  The  higher  current  density  on  the  inner  wall  of  the  bend  causes  the  inner  region  of  the  band  to  move 
faster  than  the  outer  region,  resulting  in  an  asymmetric  profile  and  eventual  band  broadening.  The  greater  distance  traveled 
by  the  species  along  the  outer  wall  further  adds  broadening  and  asymmetry  to  the  band.  As  it  travels  past  the  u-shaped 
bend  the  band  takes  the  shape  of  a  parallelogram.  In  principle,  this  shape  can  be  corrected  by  allowing  the  band  to  negotiate 
an  equal  and  opposite  turn  downstream  of  this  bend 
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The  simulations  for  this  configuration  at  the  same  field  strength  are  also  shown  in  Fig.  7  at  the  same  times  as  the 
experiment.  The  agreement  between  experiment  and  simulation  is  quite  good  indicating  that  the  diffusivity  and  mobility 
extracted  from  the  straight  channels  are  reasonably  valid  in  the  more  complex  curved  channels  as  well. 

We  have  also  conducted  similar  experimental/numerical  comparisons  for  electroosmotic  flow  to  extract  an 
electroosmotic  mobility  as  a  function  of  the  field  strength.  The  procedure  is  similar  to  the  electrophoretic  case  and  is  not 
reproduced  here. 

The  net  result  of  the  work  presented  here  is  quite  encouraging.  Using  experimental  data  generated  in  fairly  simple  test 
structures,  we  are  able  to  generate  the  parameter  required  for  the  numerical  analysis  of  more  complex  structures.  Simple 
straight  channel  structures  such  as  those  used  in  the  baseline  cases  are  generally  easy  to  manufacture,  set  up  and  test 
experimentally.  A  test  procedure  of  this  nature  can  also  be  used  to  ensure  conformity  of  manufactured  structures  to  some 
baseline  specifications  for  the  device.  This  is  analogous  to  procedures  used  in  both  the  MEMS  and  IC  areas. 

7.  CONCLUSIONS 


We  have  developed  Netflow,  a  CAD  tool  for  chemical  transport  in  microchannels.  The  tool  is  built  up  on  the  existing 
framework  for  MEMCAD,  which  is  a  general  purpose  tool  for  the  design  of  MEMS  and  microfluidic  systems.  Netflow 
enables  the  simulation  of  electrophoretic,  electroosmotic,  diffusive  and  pressure-driven  flow  in  microchannels.  In  its 
application  to  practical  design  problems,  however,  CAD  tools  require  the  knowledge  of  material  and  physical  parameters, 
that  are  generally  not  known.  In  this  paper,  we  demonstrate  a  mechanism  to  first  extract  the  parameters  of  interest  from 
simple  test  structures  and  then  apply  them  in  the  analysis  of  more  complex  structures.  This  procedure  is  analogous  to  that 
used  in  both  MEMS  and  IC  design. 

The  Simulation  Manager  tool  in  Memcad  allows  the  computation  of  the  effect  of  parametric  variation  in  a  given 
problem.  This  procedure  is  employed  in  iterating  over  a  range  of  values  of  the  parameter  of  interest,  and  then  comparing 
with  experimental  measurements  to  extract  the  required  parameter  value.  We  have  demonstrated  this  procedure  for  the 
extraction  of  diffusivity  and  electrophoretic  mobility  data  and  have  conducted  measurements  for  electroosmotic  mobility  as 
well.  Using  data  obtained  in  straight  channels  we  were  able  to  make  reasonably  good  predictions  of  flow  through  curved 
channels,  thereby  validating  this  approach.  The  next  step  is  to  implement  this  procedure  in  the  design  of  a  practical  micro¬ 
device. 
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Figure  1  :  Schematic  of  Design  Process  employed  in  Material  Property  Extraction. 
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Figure  3  :  Graph  of  the  simulated  band  width  at  two  distinct  times  (t  =  0.99 
sec,  triangles  and  t  =  1 .96  sec,  squares)  vs.  a  range  of  predicted  diffusivities. 
Circles  correspond  to  experimental  values  for  band  width  which  were  used  to 
extract  the  diffusivity. 
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Figure  4  :  Left:  Experimental  electrophoretic  flow  in  an  elastomer  channel,  E  =  85.8 
V/cm.  Right:  Column  integration  tool  from  NetFlow  used  to  form  2-D 
image/histogram  corresponding  to  the  fluorescence. 
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Fig  5.  Velocity  vs.  field  strength  used  to  obtain  the  mobility  of  uncaged 
fluorescein 


Fig  6.  Schematic  (scaled  to  size)  of  the  various  glass  serpentine 
channels 
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Figure  7:  A  series  of  timed  images  of  a  fluorescent  plug  under  electrophoretic  transport  in  a 
curved  glass  channel.  Channel  diameter  =  50  mm.  Field  strength  =  76.3  V/cm.  All  time  in 
seconds. 
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ABSTRACT 

Most  microfluidic  systems  rely  on  one  of  two  manners  of  fluid 
transport:  pressure-driven  or  electrokinetically-driven  flow.  This 
investigation  focuses  on  describing  these  flows  in  microfabricated 
channels  and  small  diameter  capillary  tubes.  Flow  characterization  is 
accomplished  by  interrogation  of  micron-scale  fluid  regions  through  a 
powerful,  non-intrusive  flow  imaging  technique.  Interesting 
phenomena  have  been  observed  from  these  detailed  examinations. 
Our  results  are  presented  in  conjunction  with  an  evaluation  of 
mechanisms  that  potentially  explain  observed  deviations  from  the 
Helmholtz-Smoluchowski  equation.  In  particular,  we  show  that 
observed  perturbations  of  electrokinetic  flow  in  open  capillaries  might 
be  caused  by  induced  pressure  gradients.  We  also  show  how  these 
induced  pressure  gradients  may  globally  perturb  the  flow  in  an 
electrokinetically-driven  microfluidic  system. 

NOMENCLATURE 


Symbol 

Definition 

Units  (SI) 

E 

electric  field 

V  m"1 

Ex 

axial  electric  field 

Vnf1 

F 

Faraday  constant  (Nae) 

C  mof1 

Io 

Modified  Bessel  function 

L 

length  of  capillary 

m 

P 

Pressure 

Pa 

Q 

Volumetric  flow  rate 

m3  s'1 

R 

universal  gas  constant 

J  mof1  K' 

S 

source  term 

Wm"3 

T 

temperature 

K 

To 

wall  temperature 

K 

U 

electroosmotic  or 
electrophoretic  velocity 

m  s'1 
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u0 

electrophoretic  velocity 
at  T0 

m  s 

Co 

ion  concentration 

mol  m"3 

dV 

differential  volume 

*rt3 

m 

k 

thermal  conductivity 

W  m"1  K- 

r 

radial  dimension 

m 

to 

capillary  inner  radius 

m 

u 

velocity  distribution 

m  s'1 

X 

axial  dimension 

m 

Xo 

axial  position  of  change  in 
wall  zeta  potential 

m 

y 

distance  from  planar  surface 

m 

z 

charge  number 

E 

permittivity 

C  W1  m' 

<t> 

electric  potential 

V 

^■D 

Debye  length  . 

m 

It 

viscosity 

Pa  s 

ho 

viscosity  of  water  at  T0 

Pa  s 

Pe 

electric  charge  density 

Cm’3 

P 

resistivity 

ohm  m 

C 

zeta  potential 

V 

Ceo 

electroosmotic  zeta  potential 

V 

^EP 

electrophoretic  zeta  potential  V 

INTRODUCTION 

Originally  emerging  from  the  field  of  analytical  chemistry, 
micro  total  analysis  systems  (pTAS)  have  recently  enjoyed 
wide  interest.  Micro  total  analysis  systems  and  other 
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microfluidic  systems  rely  on  both  pressure  and  electrokinctic 
mechanisms  for  fluid  transport.  Complicated  relationships  can 
arise  between  the  physical  characteristics  of  the  microchanncls 
and  the  behavior  of  the  multi-component  fluids  flowing  through 
the  channels.  These  relationships  have  not  been  systematically 
studied  and  are  not  yet  completely  understood.  In  practice, 
researchers  are  forced  to  rely  on  the  use  of  trial  and  error  in  the 
design  of  miniaturized  fluidic  systems. 

In  this  investigation,  we  employed  a  caged  fluorescence 
imaging  technique  in  an  effort  to  determine  the  mechanisms 
that  contribute  to  the  fluid  flow  in  microfluidic  systems.  Aided 
with  characterizations  of  relatively  simple  flows,  we  endeavor 
to  accurately  predict  the  flow  behavior  in  more  complicated 
flows  and  geometries. 


BACKGROUND 

Electrokinetic  transport  refers  to  the  combination  of 
electroosmotic  and  electrophoretic  transport.  For  a  historical 
review  of  electrokinetic  theory,  see  Burgreen  and  Nakach 
(1964).  Included  below  is  a  brief  description  of  each  type  of 
transport. 


Electroosmosis 

Electroosmosis  refers  to  the  bulk  movement  of  an  aqueous 
solution  past  a  stationary  solid  surface,  due  to  an  externally 
applied  electric  field.  Electroosmosis  requires  the  existence  of 
a  charged  double-layer  at  the  solid-liquid  interface.  This 
charged  double  layer  results  from  an  attraction  between  bound 
surface  charges  and  ions  in  the  passing  fluid.  In  glass 
capillaries,  surface  silanol  groups  become  deprotonated  and. 
therefore,  are  negatively  charged.  This  negatively  charged 
surface  attracts  positive  ions  present  in  the  flow.  In  the 
situations  addressed  in  this  paper,  only  a  very  thin  layer  near 
the  wall  has  a  net  charge. 

Rice  and  Whitehead  (1965)  give  a  complete  analysis  of 
electroosmotic  flow  in  round  capillaries;  here  we  quote  only  a 
few  significant  results  to  which  we  will  refer  later.  At 
equilibrium,  a  Boltzmann  distribution  of  the  charges  in  the 
solution  will  exist  (Probstein,  1994)  with  the  charge  density 
distribution  represented  by: 
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The  potential  field  is  determined  by  solving  the  Poisson- 
Boltzmann  equation,  resulting  in  the  following  charge  density, 
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In  our  experiments,  Ttyr,,  is  approximately  10'4;  thus  the  double 
layer  is  very  thin.  The  equation  of  motion  for  steady,  low 
Reynolds  number,  incompressible  flow  is: 


//  V2U  =  -yC>£E  +  VP  .  (3) 


In  the  limit  of  a  small  Debye  length,  A,D,  solving  Eq.  (3)  yields 
the  Helmholtz-Smoluchowski  equation  for  the  electroosmotic 
velocity. 
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The  Helmholtz-Smoluchowski  equation  predicts  a  “plug- 
like"  velocity  profile  when  the  Debye  length  is  much  less  than 
the  capillary  diameter.  The  electroosmotic  velocity,  U,  is 
approximately  0.1  mm/sec,  for  £  =  0.1  V,  Ex  =  100  V/cm,  and  p 
equal  to  the  viscosity  of  water,  (Probstein,  1994). 


Electrophoresis 

Electrophoresis  describes  the  motion  of  a  charged  surface 
submerged  in  a  fluid  under  the  action  of  an  applied  electric 
field.  Considering  the  case  of  a  charged  dye  molecule,  it  can  be 
shown  (Probstein,  1994)  that  the  electrophoretic  velocity  of  the 
dye  is  again  described  by  the  Helmholtz-Smoluchowski 
equation. 


T  &  (j  n -p  Er 

jj  _  _b_EP - x_  (5) 

P 

The  electroosmotic  zeta  potential  is  a  property  of  the  capillary 
surface  while  the  electrophoretic  zeta  potential  is  a  property  of 
the  charged  dye.  In  general,  these  two  zeta  potentials  will  not 
be  the  same. 

It  is  important  to  note  that  electroosmosis  results  in  a  net 
mass  transfer  of  the  aqueous  solution;  whereas,  electrophoresis 
causes  movement  of  charged  particles  or  molecules  through  a 
stationary  solution.  For  a  net-neutral  solution  of  charged 
molecules,  electrophoresis  will  not  result  in  bulk  motion  of  the 
solvent.  Also  note  that  both  derivations  shown  above  result  in 
a  velocity  that  is  not  a  function  of  the  radial  position;  hence,  it 
is  termed  a  “plug-like"  flow. 


EXPERIMENTAL  DETAILS 

Two  different  types  of  microchannels  were  used  in  this 
investigation.  Initially,  we  used  rectangular  silicone 
microchannels.  We  then  began  investigating  flows  within 
round  glass  capillaries  with  an  internal  diameter  of  100  pm. 
The  capillaries  were  employed  in  an  effort  to  isolate 
unexplained  flow  phenomena  observed  in  our  initial  silicone 
microchannel  investigation. 

To  create  the  silicone  channels,  silicon  microchannel 
molds  were  fabricated  as  “negative”  masters  using  standard 
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photolithography.  An  STS  deep  reactive  ion  etch  (DRIE)  was 
used  to  create  raised  rectangular  structures  (inverted  channels) 
with  dimensions  of  10-200  pm  in  width,  5  cm  in  length  and  40 
pm  in  height.  An  RTV  silicone  (Dow  Coming  Sylgard  184) 
was  used  to  form  positive  replicas  of  the  silicon  structures 
(Effenhauser,  et  al.  1997).  A  photograph  of  a  sealed  replica  is 
shown  in  Fig.  1.  The  elastomer  channels  were  sealed  with 
either  glass  microscope  slides  or  elastomer  coated  glass 
microscope  slides. 


Figure  1.  Photograph  of  the  elastomer  device  with  three 
test  microchannels.  The  circles  are  the  injection  wells  and 
the  squares  are  used  to  study  pressure-driven  flow  in  a 
sudden  expansion. 

Our  experimental  approach  utilized  caged  fluorescent  dyes 
to  image  the  flows.  These  dyes  were  activated  (uncaged)  by  a 
UV  laser  beam  and  tracked  by  fluorescence  imaging,  see  Fig.  2, 
(Paul  et  al,  1998).  In  contrast  to  standard  methods  injecting 
dye,  this  visualization  technique  allows  definition  of  narrow, 
fluorescent  regions  at  virtually  any  location  along  the  channel. 
Specifically,  this  allows  for  a  precise  definition  of  a  small 
fluorescent  region  and  therefore  an  improved  determination  of 
the  dye  profile.  Therefore,  this  caged  fluorescence  technique 
has  advantages  over  injection  methods  that  must  transport  the 
dye  from  the  location  of  injection  to  the  region  of  interest. 

The  caged  fluorophores  were  purchased  from  Molecular 
Probes  and  used  in  millimolar  quantities.  For  the  pressure- 
driven  studies,  a  caged  fluorescein  was  dissolved  in  distilled 


UV  pulse  initiates  dye  line 


Figure  2.  Top:  Schematic  of  the  experiment.  The  Nd:YAG 
laser  is  used  to  uncage  the  dye  and  a  Microblue  473  nm 
laser  is  used  to  excite  the  fluorescence.  Bottom: 
Illustration  of  the  uncaging  process.  The  uncaging  laser 
defines  the  starting  fluorescent  fluid  volume. 

water.  For  the  electrokinetic  studies,  a  caged  rhodamine  was 
dissolved  in  Tris-EDTA  buffer.  The  dyes  are.  negatively 
charged  when  uncaged  and  may  be  either  charged  or  uncharged 
when  caged.  When  caged,  the  dye  solution  is  non-fluorescent. 
Upon  photo-activation  of  a  UV  light  illuminated  volume,  the 
protecting  group  is  cleaved  and  the  dye  becomes  fluorescent. 
Approximately  100  pJ  (per  pulse)  of  355  nm  light  from  a 
pulsed  Nd:YAG  laser  uncages  the  dye.  The  uncaged  region  is 
roughly  20  pm  along  the  length  of  the  channel  or  capillary  and 
extends  across  its  full  width  (Fig.  2).  The  laser  beam  was 
focused  using  a  UV  microscope  objective  that  defined  a  sharp 
start  zone.  A  continuous  wave  Microblue  diode  pumped  laser  at 
473  nm  (Uniphase)  was  used  for  fluorescence  excitation  after 
activation  (Fig.  2).  Fluorescence  images  of  the  molecules 
diffusing  and  moving  in  the  local  flow  were  collected  using  a 
microscope  objective  and  a  video  rate,  interlaced  camera 
(Texas  Instruments  MC-780PH). 
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RESULTS 

We  have  investigated  transport  in  microfluidic  systems  and 
present  experimental  studies  of  pressure-driven  and 
electrokinetic  flows.  We  also  present  numerical  simulations 
which  were  performed  using  the  NetFlow  module  in  Memcad 
(MEMCAD,  1998).  The  NetFlow  module  employs  a  three- 
dimensional  finite  element  based  tool  to  solve  the  Navier- 
Stokes  equations.  We  note  that  the  pressure-driven  flow  data 
show  good  agreement  between  the  preliminary  experiments  and 
the  simulations.  Measurements  of  the  electrokinetic  flows, 
however,  reveal  deviations  from  simple  theory. 

The  evolution  of  a  pressure-driven  fluorescent  profile  is 
compared  to  a  corresponding  sequence  of  simulations  in  Fig.  3. 
The  agreement  is  good,  but  note  that  there  is  dye  retention  near 
the  walls  that  is  not  duplicated  in  the  numerical  simulations. 
This  dye  retention  could  be  due  to  surface  roughness.  Atomic 
force  microscopy  indicates  that  the  silicone  channel  has  an 
RMS  surface  roughness  of  ~  7  nm. 


Figure  3.  Experimental  (left)  and  simulated  (right)  pressure- 
driven  flow  entering  a  sudden  expansion  in  an  elastomer 
structure.  The  channel  entering  the  expansion  is 
rectangular  (100  pm  wide  and  40  pm  deep) 


The  evolution  of  an  electrokinetically-driven  fluorescent 
profile  in  a  round,  fused-silica  capillary  is  shown  in  Fig.  4. 
These  images  reveal  that  the  velocity  profile  appears  to  be 
constant  across  the  channel,  in  agreement  with  Eq.  (4).  Figure 
5  shows  electrokinetic  flow  in  a  DBwax-coated  (polyethylene 
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Figure  4.  Electrokinetic  flow  in  a  75  pm  diameter  fused- 
silica  capillary  exhibiting  a  plug-like  velocity  profile. 
Applied  field  strength  was  350  volts/cm.  The  numbers  with 
each  image  correspond  to  the  time  in  milliseconds  after  the 
uncaging  pulse. 
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Figure  5.  Electrokinetic  flow  in  DBwax-coated  (polyethylene 
glycol),  100  pm  diameter  fused-silica  capillary.  The  wax- 
coated  capillaries  were  purchased  from  J&W  Scientific. 
Applied  field  strength  was  245  volts/cm.  The  numbers  with 
each  image  correspond  to  the  time  in  milliseconds  after  the 
uncaging  pulse. 


glycol)  capillary.  The  wax  coating  is  used  to  suppress 
electroosmotic  flow,  so  that  the  movement  of  charged  particles 
in  the  flow  is  due  to  electrophoresis  only.  However,  these 
images  show  a  pronounced  parabolic  velocity  in  the  direction 
opposing  the  electric  field.  This  profile  does  not  agree  with  Eq. 
(5),  and  would  be  difficult  to  image  without  a  high-resolution 
measurement  technique  such  as  caged  florescence  imaging. 
Note  that  we  observed  no  motion  of  the  dye  without  application 
of  the  electric  field,  verifying  that  there  was  no  externally 
applied  pressure  gradient.  Parabolic  components  to  the 
electrokinetic  velocity  profile  have  previously  been  observed  in 
capillaries  (Paul  et  al,  1998)  and  in  microchannels  (St.  John  et 
al,  1998). 
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DISCUSSION 

In  Fig.  4,  we  see  that  the  imaged  velocity  profiles  for 
electrokinetic  flow  are  in  good  agreement  with  Eq.  (4).  The 
flow  is  a  combination  of  electroosmosis  and  electrophoresis, 
and  the  theory  outlined  above  shows  that  both  components 
yield  a  flat  profile.  Figure  5,  however,  shows  significant 
deviations  from  a  flat  profile.  This  flow  suggests  a  velocity 
profile  similar  to  that  observed  in  a  pressure-driven  flow. 

In  the  following  sections,  we  explore  two  mechanisms  that 
potentially  explain  this  observed  departure  from  plug-like  flow: 
variations  in  viscosity  due  to  temperature  gradients,  and 
induced  pressure  gradients  caused  by  a  non-uniform  wall 
potential  (i.e.  non-uniform  wax  coating  in  the  case  of  Fig.  5). 


u(T)  =  i£m dkm  (9) 

M(T) 

The  predicted  velocity  profiles  for  both  an  open  capillary 
and  a  DBwax-coated  capillary  each  show  an  expected  radial 
velocity  variation  of  less  than  0.2%  (Fig.  6).  The  radial 
velocity  distribution  is  parabolic  in  nature,  but  its  magnitude  is 
too  small  to  explain  the  dye  profiles  seen  in  Fig.  5.  To  obtain  a 
significant  radial  variation  in  the  electrophoretic  velocity,  one 
would  need  to  increase  the  field  strength  or  the  buffer 
concentration  considerably.  Our  temperature  variation  model 
and,  hence  these  conclusions,  are  in  close  agreement  with  those 
presented  by  Grushka,  et  al.  (1989). 


Variation  of  viscosity 

In  this  section,  we  examine  the  possibility  that  heat 
generation  within  the  capillary  could  cause  the  electrophoretic 
velocity  profile  to  deviate  from  plug-like.  To  this  end,  we 
investigate  the  temperature,  viscosity  and  corresponding 
velocity  distributions  across  the  capillary. 

In  cylindrical  coordinates,  we  write  the  temperature 
governing  equation  for  our  system  as: 


ld_ 
r  dr 
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The  source  term,  S,  has  units  of  power  per  unit  volume.  In  this 
analysis,  we  assume  p  remains  a  constant  (or  at  least  not  a 
strong  function  of  temperature)  based  on  the  experimentally 
determined  fluid-filled  capillary  resistance  and  capillary 
dimensions.  This  empirically  determined  resistivity  agrees 
closely  with  those  tabulated  in  the  literature  (Masliyah  1994). 
Integrating  twice  and  applying  the  boundary  conditions  of  finite 
temperature  at  the  centerline  and  a  constant  wall  temperature  of 
T0,  we  obtain  the  temperature  distribution 


T(r)  = 


p4k 
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As  an  initial  approximation,  we  assume  that  the  fluid  within  our 
capillary  has  the  same  thermal  properties  and  viscosity  as 
water.  As  such,  we  have  chosen  to  model  the  viscosity  as  a 
function  of  temperature, 
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with  T0  =  293. 15K  and  the  constants  given  by  a  =  -1.94,  b  =  - 
4.80,  and  c  =  6.74,  as  suggested  by  Reid  (1987). 

In  order  to  relate  the  temperature-dependent  viscosity  to 
the  electrophoretic  velocity  profile  within  the  capillary,  we 
return  to  our  analysis  of  electrophoresis.  Utilizing  Eq.  (5)  and 
neglecting  the  temperature  dependence  of  e,  or  Ex,  we  arrive 
at: 
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Figure  6.  Electrophoretic  velocity  distribution  for  an  open 
capillary  and  a  DBwax-coated  capillary 


From  this  analysis,  we  conclude  that  the  viscosity  and, 
therefore,  velocity  do  not  vary  significantly  across  the  capillary 
for  the  prescribed  small-bore  capillary  parameters,  buffer 
concentrations,  and  field  strengths.  Therefore,  Joule  heating 
induced  velocity  variation  does  not  appear  to  be  a  feasible 
mechanism  for  the  observed  deviations  from  plug-like  flow. 


Induced  Pressure  Gradients 

In  deriving  Eq.(4),  the  pressure  gradient  term  was 
neglected  because  there  were  no  externally  applied  pressure 
gradients.  However,  this  expression  will  not  be  correct  if  the 
electrokinetic  transport  induces  a  pressure  gradient.  We  now 
explore  how  the  electrokinetic  flow  could  potentially  induce 
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pressure  gradients.  If  we  take  the  divergence  of  Eq.  (3),  and 
apply  continuity  for  the  fluid,  (V  ■  u )  =  0 ,  we  derive  the 
following  expression. 

W2P  =  E-VpE +(V -E)pe  (10) 

Equation  (10)  demonstrates  how  the  pressure  field  evolves 
from  gradients  in  the  charge  density  and  applied  field. 
Referring  to  Eq.  (2),  we  see  that  the  first  term  on  the  right  hand 
side  of  Eq.  (10)  will  be  non-zero  if  the  wall  potential  changes  in 
the  axial  direction.  Nonuniformities  in  the  wall  potential  could 
exist  due  to  contamination  in  the  capillary,  variations  in  wall 
coatings,  or  gradients  in  the  buffer  pH.  Using  classical 
techniques  for  solving  Poisson’s  equation  (e.g.  Green’s 
functions),  it  should  be  possible  to  develop  a  general  solution  to 
Eq.  (10)  for  an  arbitrary  zeta  potential  distribution  along  the 
wall  of  a  cylindrical  capillary.  In  this  paper,  a  direct  method 
will  be  used  to  solve  for  the  pressure  field  in  a  capillary  with  a 
wall  potential  as  shown  in  Fig.  7. 
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capillary  is  a  known  reference  value  (i.e.  atmospheric  pressure). 
For  an  incompressible  liquid,  the  volumetric  flow  rates  must  be 
the  same  in  regions  1  and  2.  Furthermore,  we  know  that  the 
pressure  is  continuous  and  piecewise  linear  since  we  assumed 
that  the  velocity  profile  does  not  vary  axially  in  either  region 
(i.e.  the  pressure  gradient  must  be  constant  in  both  regions). 
Therefore,  we  arrive  at  the  following  expression  for  the 
pressure  in  the  capillary',  where  x  is  the  axial  position  and  x0  is 
the  axial  location  of  the  discontinuity  in  wall  potential. 
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Figure  7.  Schematic  of  capillary  with  step  change  in  wall 
potential.  The  anticipated  velocity  profiles  are  also  shown. 


Figure  8.  Non-dimensional  pressure  in  a  capillary  with  a 
step  change  in  wall  potential. 


Let  us  assume  that  the  fluid  velocity  is  only  a  function  of 
radial  position  in  both  Region  1  and  Region  2  shown  in  Fig.  7. 
Note  that  Region  1  describes  a  section  of  the  capillary  with  zero 
wall  potential;  whereas,  Region  2  has  a  finite  wall  potential.  In 
other  words,  we  will  assume  that  axial  velocity  gradients  are 
confined  to  a  vanishingly  small  region  near  the  discontinuity  in 
wall  potential.  With  this  assumption,  the  volumetric  flow  rates 
in  regions  one  and  two  are  shown  below. 
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We  consider  the  case  where  the  pressure  at  either  end  of  the 


P(x)  = 


*ExeC 
..  2  1 


j{L-x„),  x<xo 
~{L~x),x>x0 


(12) 


The  pressure  field  corresponding  to  Fig.  7  is  plotted  in  Fig. 
8.  We  see  that  the  pressure  gradient  is  constant  in  each  section 
and  acts  to  ensure  that  the  mass  flow  rates  are  equal  in  each 
region.  Figure  9  shows  the  predicted  velocity  profiles  in 
regions  1  and  2.  The  velocity  in  the  first  section  is  due  to 
pressure  alone,  since  the  wall  potential  there  is  zero.  The 
velocity  profile  in  the  second  section  is  simply  a  superposition 
of  the  flat  electroosmotic  profile  and  the  Poiseuille  profile  due 
to  the  induced  pressure  gradient.  The  pressure  gradient  in 
each  region  is  readily  calculated  from  Eq.  (12).  The  velocity 
profiles  in  Fig.  9  were  also  shown  schematically  in  Fig.  7. 
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Figure  9.  Radial  velocity  profiles  in  region  1  and  region  2, 

To  understand  the  solution  qualitatively,  one  can  imagine 
that  the  electroosmotic  flow  in  Region  2  “pulls”  the  flow  in 
Region  1.  In  other  words,  the  electroosmotic  movement  of 
liquid  in  Region  2  tends  to  create  suction  at  the  interface 
between  the  two  regions.  Hence,  the  pressure  drops  at  the 
interface  until  the  mass  flow  rates  in  each  region  are  equal. 
Note  that  the  pressure-driven  component  subtracts  from  the 
electroosmotic  flow  in  Region  2  while  only  pressure-driven 
flow  is  present  in  Region  1.  The  analytically  calculated 
velocity  profiles  also  agree  qualitatively  with  the  numerical 
simulation  shown  in  Fig.  10  (MEMCAD,  1998).  The 
simulation  confirms  that  axial  gradients  in  the  velocity  field  are 
confined  to  a  region  less  than  a  single  diameter  from  the  wall 
potential  discontinuity. 

Superpositions  of  Eq.  (12)  can  be  used  to  find  the  pressure 
field  caused  by  one  or  many  small  regions  of  zero  wall 
potential.  This  is  then  a  model  of  a  “dirty”  silica  capillary. 
Figure  1 1  demonstrates  the  nondimensional  pressure  field  in  a 
capillaiy  with  three  “dirty”  regions,  each  covering  1%  of  the 
capillary’s  length. 

Note  that  the  pressure  gradients  in  the  clean  regions  are 
identical  and  that  the  localized  discontinuities  in  wall  potential 
perturb  the  flow  everywhere.  Equation  (12)  can  also  be 
superposed  to  find  the  pressure  field  in  a  coated  capillary  (zero 
wall  potential)  with  small  “clean”  regions  where  the  wall 
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Figure  10.  Simulation  of  flow  in  a  capillary  with  a  step 
change  in  wall  potential. 
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Figure  11.  Nondimensional  pressure  field  in  a  capillary 
with  three  “dirty”  spots 

potential  is  non-zero.  Figure  12  shows  the  resulting 
velocity  profile  in  the  coated  regions,  for  various  amounts  of 
exposed  (uncoated)  capillary.  The  induced  pressure  gradient 
model  presents  a  rational  explanation  of  the  flow  imaged  in  Fig. 
5,  as  it  admits  a  parabolic  profile  for  reasonable  variations  in 
wall  potential. 
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Figure  12.  Velocity  profile  in  a  coated  region  of  a  wax 
coated  capillary  with  small  regions  uncoated.  The  wax 
coating  acts  to  prevent  electroosmotic  flow. 
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CONCLUSIONS 

We  have  imaged  pressure-driven  and  electrokineticallv- 
driven  flows  in  silicone  microchannels  and  silica  capillaries 
using  a  powerful  caged  fluorescence  technique.  The  pressure- 
driven  flows  agree  with  theory  and  numerical  simulations. 
Simple  theory  predicts  that  the  electrokinetic  flows  should 
display  uniform  (plug-like)  velocity  profiles.  However,  in 
some  cases,  we  observe  large  parabolic-like  components  to  the 
electrokinetic  velocity  profiles. 

We  explored  two  mechanisms  in  an  attempt  to  explain  the 
observed  velocity  profiles:  viscosity  variations  and  induced 
pressure  gradients.  We  have  demonstrated  that  Joule  heating 
and  the  subsequent  temperature  and  viscosity  gradients  do  not 
adequately  explain  the  observed  parabolic  velocity  profile,  but 
that  induced  pressure  gradients  are  a  possible  explanation.  An 
important  consequence  of  the  induced  pressure  gradient 
analysis  is  that  small  regions  of  non-uniform  wall  potential 
create  pressure  gradients  everywhere  in  the  flow.  This 
conclusion  has  important  implications  for  those  wishing  to 
build  lab-on-a-chip  devices,  since  pressure  gradients  tend  to 
increase  the  dispersion  of  sample  packets  in  the  flow.  We  are 
currently  performing  experiments  that  investigate  the  proposed 
induced  pressure  gradient  mechanism. 
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Abstract 

Design  analysis  of  PCR  (Polymerase  Chain  Reaction) 
microwells  are  conducted  to  observe  the  contamination 
between  neighboring  wells  on  a  PCR  chip  during  the 
reaction.  Thermal  effects  are  incorporated  to  represent  the 
temperature  cycling  characteristic  of  PCR.  The  analyses  are 
intended  to  be  predictive  -  and  to  be  applied  in  the  design  of 
high  throughput  PCR  devices.  Parametric  analyses  were 
conducted  to  observe  the  effect  of  the  allowable  geometry 
variations  on  the  contamination  of  the  neighboring  wells, 
and  the  effect  of  the  variation  in  the  diffUsivities  of  the 
constituent  chemical  components.  The  analyses  help  in 
developing  a  set  of  curves  that  can  be  used  to  develop  a 
model  to  represent  the  behavior  of  the  PCR  well.  This 
model  can  then  be  implemented  in  system  models  for  the 
entire  device. 

Introduction 

There  is  a  wide  interest  in  micron-scale  integrated 
chemical/biochemical  analysis  or  synthesis  systems,  also 
referred  to  as  lab-on-a-chip  [1-3].  On-chip  analysis  systems 
consist  of  active  and  passive  components.  Examples  of  the 
former  are  reaction  chambers,  switching  joints  etc.,  whereas 
the  latter  include  transport  and  separation  channels.  The 
behavior  of  these  components  is  generally  dependent  on  a 
number  of  parameters,  including  flow  conditions  and  fluid 
properties.  Researchers  are  forced  to  use  costly  trial-and- 
error  experimental  methods  to  understand  and  design  such 
microfluidic  systems. 

CAD  analyses,  on  the  other  hand,  can  be  a  valuable  aid 
in  the  design  of  microfluidic  systems.  Numerical  analyses 
provide  significant  insight  into  the  fluid  mechanics  in  these 
systems.  They  also  allow  the  extraction  of  material  and  flow 
properties  that  are  generally  not  well  documented.  Such 
tools  help  the  designer  to  explore  a  much  larger  space  of 
designs  than  is  easily  available  from  experiment,  and  do  so 
in  a  quantitative  way  which  enables  the  extraction  of  key 
parameters  for  improved  or  optimal  operation  of  common 
microchemical  system  components.  Finally  numerical 
analyses  enable  the  extraction  of  component  models  for 
each  individual  component  in  a  system  model  of  the  device. 


Several  researchers  have  reported  CAD-based  analyses 
of  microfluidic  components  [4-6],  These  include 
components  used  in  injection,  transport,  as  well  as  mixing 
and  reactions.  These  authors  [7-8]  have  also  reported 
coupled  experimental  and  numerical  analyses  for  various 
microfluidic  components.  The  component  models  can  be 
employed  in  understanding  the  behavior  of  the  entire  device 
as  a  system.  The  value  of  the  system  model  in 
understanding  the  performance  of  the  system  is  dependent 
on  the  ability  of  the  component  models  in  representing  the 
important  physics  in  the  component. 

An  example  of  a  microfluidic  device  is  shown  in  Fig.  1 . 
The  device  is  designed  for  high  throughput  polymerase 
chain  reactions  (|iPCR).  The  primary  components  in  this 
device  are  the  microwells  and  the  interconnects  between 
neighboring  wells  and  to  the  device  loading  site.  The 
devices  are  operated  by  cycling  the  temperature  to  control 
the  PCR  reaction.  The  behavior  of  the  device  is 
characterized  by  the  ability  of  the  wells  to  contain  the 
sample  during  the  reaction  with  no  leakage  into  the 
neighboring  wells.  The  containment  is  affected  significantly 
by  the  geometry  of  the  interconnects  between  the  wells  and 
the  effect  of  the  temperature  cycling  on  the  diffusivity  of  the 
sample.  Good  design  of  the  interconnects  aims  to  minimize 
the  contamination  of  neighboring  wells  during  the  reaction. 


Figure  1:  Geometry  of  Current  PCR  Chip 
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Figure  2:  PCR  Microwell  Configurations  - 
Diagonal  (left)  and  Side-by-side  (right) 


In  general,  the  entire  space  of  design  parameters  is  fairly 
difficult  to  analyze  experimentally.  Numerical  methods  on 
the  other  hand  are  a  viable  mechanism  -  and  are  the  focus  of 
this  paper. 

We  present  here  a  numerical  analysis  of  the  PCR  wells 
during  the  reaction.  The  effect  of  the  temperature  cycling  on 
the  diffusivity  is  considered  in  the  analysis.  The  effect  of 
varying  the  cross  sectional  area  of  the  interconnects  is  also 
presented  to  demonstrate  the  potential  of  numerical  methods 
to  effectively  analyze  the  design  space  for  these  problems. 
Finally,  we  will  discuss  a  mechanism  of  developing  a 
system  model  for  the  device. 

Device  Geometry 

The  geometry  of  the  current  plastic  PCR  chip  is  shown  in 
Figure  1.  It  consists  of  an  array  of  49,  1 -microliter  wells.  The 
nearest  neighbor  periodicity  in  the  well  layout  can  have  two 
configurations  on  the  chip,  as  shown  in  Figure  2:  side-by-side 
(left)  or  across  the  diagonal  (right).  A  possibility  for  a 
higher  density  design  may  consist  of  an  array  of  1024  100 
nanoliter  wells  fabricated  from  silicon  using 
photolithography  and  an  anodically  bonded  pyrex  glass  cover 
to  seal  the  wells.  As  in  the  chip  in  Figure  1,  the  wells  need  to 
be  connected  to  a  microchannel  running  the  length  of  the 
chip  and  column  and  connecting,  via  a  cross  channel,  to  an 
inlet  port. 

The  PCR  chip  shown  in  figure  1  is  loaded  with  different 
reagents  by  drying  them  down  in  each  well  using  a 
dispensing  robot  [9].  The  back  is  then  sealed  using  a 
polymer  tape  and  the  wells  are  evacuated  and  then  loaded 
with  sample  via  the  inlet  port.  The  device  is  temperature 
cycled  on  a  modified  research  instrument  built  for  real  time 
detection  of  the  (iPCR  reaction.  Details  of  the  temperature 
cycling  are  shown  Fig.  3.  At  the  outset  the  system  is  at  room 
temperature  (30  °C).  The  temperature  is  then  raised 
uniformly  at  4  °C/sec  from  room  temperature  to  95  °C  and 
held  for  1 0  seconds.  It  is  then  lowered  at  the  same  rate  to  60 


Time  (sec) 

Fig.  3:  Temperature  Cycling  in  ftPCR  wells. 

°C  and  held  for  5  seconds.  The  process  then  repeats  for  a 
total  reaction  time  of  20  minutes.  Note  that  the  cycling 
shown  in  Fig.  3  repeats  for  the  specified  time.  Cycling  times, 
in  general,  are  specific  for  a  given  device  material,  volume  of 
sample,  and  reagent.,  The  times  used  here  are  characteristic 
for  microwells  and  were  used  in  the  simulation. 

Numerical  Analysis 

Predictive  CAD  analyses  for  the  PCR  well 
configurations  were  conducted  to  observe  the  effect  of  the 
temperature  cycling  on  the  containment  of  the  sample 
within  the  well.  The  leakage  into  the  nearest  neighbor  well 
was  also  analyzed  in  the  process.  The  different  allowable 
configurations  are  also  studied  as  part  of  the  analysis.  The 
numerical  analyses  were  conducted  using  FlumeCAD. 
FlumeCAD  is  an  integrated  design  environment  consisting 
of  3D  design,  modeling  and  simulation  software  tools  which 
enable  the  creation  of  complex  microfluidic  devices. 
Netflow,  which  is  a  part  of  the  suite,  is  specific  to  the 
analysis  of  chemical  transport  in  microchannels.  It  is  an 
ongoing  program,  and  is  a  result  of  a  collaboration  between 
Microcosm,  Perkin-Elmer  Corp.  and  Stanford  University. 

The  numerical  solution  of  the  above  equations  is 
performed  using  a  three-dimensional  finite  element  based 
CFD  engine.  The  Stokes-Einstein  formulation  was  used  to 
determine  the  diffusivity  variation  with  temperature.  Over 
the  cycle  used  here  the  diffusivity  varies  by  a  factor  of  about 
3.  The  FlumeCAD  interface  allows  the  user  to  define  the 
geometry  from  layout  and  process,  or,  from  a  model  built  in 
the  solid  modeler  module.  Once  meshed,  the  model  is 
available  to  a  suite  of  solvers  for  various  physical  domains, 
and  to  the  simulation  manager  to  parametrize  the 
dependency  in  any  domain. 
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Figure  4:  Comparison  of  Experiment  and  Simulation 
for  Integrated  Mass  Leakage  out  of  region. 


Fig  5:  Solid  Model  of  PCR  Wells. 


Results 

Results  for  the  PCR  well  analysis  are  discussed  in  this 
section.  Initially  validation  of  the  numerical  model  is 
presented,  followed  by  analysis  of  the  different  well 
configurations.  Simulations  were  performed  both  at  ambient 
temperatures  as  well  as  with  temperature  cycling  to  compare 
the  difference  in  containment  with  temperature  cycling  in 
the  analysis. 

A  range  of  diffusivities  is  present  in  the  experiment, 
corresponding  to  both  reactants  and  products.  The  lowest 
diffusivity  is  5  e-7  cm2/s,  whereas  the  highest  is  1  e-5  cm2/s 
corresponding  to  the  short  chain  probe.  In  the  simulations 
we  analyzed  diffusivities  over  this  range  to  ascertain  the 
effect  of  the  various  components  on  the  containment  in  the 
microwells. 


Experimental  Validation:  To  validate  our  methods, 
diffusivity  experiments  were  conducted  in  microchannels  to 
observe  the  leakage  of  dye  out  of  a  narrow  region.  The 
experiments  were  conducted  using  fluorescein  dye  in  a 
straight  section  of  a  microcharmel  (not  shown).  The 
comparison  between  the  experiment  and  simulation,  is 
shown  in  Fig.  4,  and  was  conducted  by  digitizing  the  images 
and  integrating  the  normalized  intensity  over  a  specific 
region  -  the  integrated  intensity  corresponds  to  the  amount 
of  dye  in  the  specific  region.  The  agreement  between 
experiment  and  simulation  is  good  over  the  entire  time 
sequence,  as  shown  in  the  figure.  This  formulation  was  used 
to  predict  the  extent  of  containment  during  PCR. 

PCR  Well  Simulations :  The  PCR  containment  analysis  is 
conducted  using  the  dimensions  from  the  next  generation 
designs.  For  the  purposes  of  the  simulation,  detailed 
analyses  were  conducted  on  a  reduced  configuration, 


Fig.  6:  Species  Concentration  of  Dye  in  Containment  Wells  over  time.  Time  increments  are  in  minutes. 
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Figure  7:  Containment  in  well  during  reaction 
and  contamination  of  Neighboring  well. 
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Figure  9:  Contamination  time  as  a  function 
of  interconnect  side  length. 
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Fig  8:  Leakage  rate  from  PCR  well  with  and  without 
temperature  cycling 


consisting  of  two  nearest  neighbor  wells.  The  solid  model 
for  this  configuration  is  shown  in  Fig.  5.  This  model  is 
adequate  to  analyze  the  effects  of  the  temperature  cycling 
and  the  diffusivity  on  the  containment/contamination 
characteristics  of  the  design  geometry.  The  problem  also 
becomes  considerably  more  tractable. 

Figure  6  shows  a  simulated  timed  series  of  diffusion  of 
material  out  of  a  microwell  and  along  an  interconnect.  The 
region  inside  the  well  is  detected  during  the  PCR  reaction 
and  therefore,  any  material  which  diffuses  out  of  the  well  in 
Figure  6  is  insignificant  providing  the  material  does  not 
diffuse  into  the  neighboring  well.  Initially,  the  reactant 
travels  out  of  the  well  quickly  as  it  establishes  a 
concentration  gradient  along  the  channel.  At  later  times,  the 
solution  moves  much  slower  as  it  diffuses  down  the  channel. 
Time  sequences  in  the  figure  are  shown  for  a  typical  PCR 
experiment  -  when  allowed  to  run  longer,  the  species 


concentrations  were  predicted  to  equalize  in  the  two  wells  in 
about  8  hours  (not  shown).  The  species  concentration 
contained  in  the  well  over  time  is  shown  in  Figure  7.  The 
concentration  entering  the  neighboring  well  is  also  shown  in 
the  figure.  These  curves  were  computed  using  a  Region-of- 
Interest  (ROI)  integration  over  the  wells  after  the 
computation.  The  ROI  integral  is  then  analogous  to  the 
collected  signal  in  a  detection  system  for  the  well.  A  curve 
such  as  Fig.  7  is  an  important  asset  in  the  design  of  these 
devices  since  it  indicates  the  allowable  reaction  times  that 
prevent  contamination  of  reagents  (and  therefore  the 
reaction)  between  neighboring  wells. 

The  effect  of  incorporating  the  temperature  cycling  is  shown 
in  Figure  8.  The  leakage  rate  from  the  chamber  is  compared 
with  the  results  of  the  simulation  at  constant  temperature. 
The  diffusion  constant  in  the  isothermal  case  was  taken  to  be 
constant  at  the  room  temperature  value.  A  different  constant 
temperature  value  would  yield  a  shift  in  the  leakage  rate,  but 
would  still  be  inaccurate  in  comparison  to  the  variable 
temperature  results.  As  the  figure  shows  the  temperature 
cycling  and  the  correspondingly  higher  diffusivity  has  a 
dramatic  effect  in  increasing  the  leakage  rate  from  the 
reaction  chamber.  This  figure  also  serves  to  highlight  the 
importance  of  including  the  relevant  physical  processes  into 
the  computation  since  the  presence  of  the  temperature  cycle 
significantly  alters  the  design  analysis  for  the  problem. 

Effect  of  interconnect  seometrv:  The  cross  sectional  area  of 
the  interconnect  can  also  have  a  significant  effect  on  the 
contamination  time  between  the  wells.  Computations  were 
conducted  using  the  Simulation  Manager  in  FlumeCAD  to 
analyze  the  effect  of  the  allowable  variations  in  the  design. 
Results  are  presented  in  Fig.  9,  showing  the  contamination 
time  for  the  different  configurations.  The  contamination  time 
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here  was  defined  as  the  amount  of  time  required  for  5  percent 
of  the  species  in  the  first  well  to  enter  the  neighboring  well. 
As  the  figure  shows,  for  the  highest  (150  pm)  case, 
contamination  occurs  well  within  the  specified  reaction  time. 
For  a  side  length  of  100  pm  as  well,  contamination  occurs  at 
approximately  10.5  minutes,  which  is  within  the  reaction 
time.  For  a  70  pm  side  length,  however,  the  contamination 
time  increases  to  over  22  minutes,  implying  that  it  is  a  more 
appropriate  choice  for  the  interconnect  cross  sectional 
dimension. 


Figure  10:  Contamination  time  as  a  function  of 
diffusivity. 


Effect  of  Diffusivity.  In  the  course  of  the  PCR  reaction 
the  diffusivities  observed  for  the  sample,  enzyme  and 
product  mixture  range  from  about  5e-7  cm2/s  to  about 
le-5  cm2/s.  The  rate  of  contamination  is  strongly 
dependent  on  the  diffusivity.  However  a  full  reaction 
analysis  for  PCR  is  computationally  intensive  -  and  is 
not  entirely  required  for  design  analysis.  We  demonstrate 
this  here,  by  studying  the  effect  of  the  diffusivity  on  the 
contamination.  In  general  the  contamination  times  lie 
between  that  for  the  highest  and  the  lowest  diffusivity  - 
the  design  can  then  be  tailored  for  an  appropriate 
specification.  The  contamination  time  has  a  linear 
dependence  on  the  diffusivity  as  seen  in  Figure  10  - 
appropriate  leakage  rates  for  specific  diffusivities  can 
then  be  extracted  from  this  curve. 


Curves  such  as  Figures  7-10  describe  the  behavior  of  the 
PCR  wells  as  a  component  in  the  overall  PCR  system.  This 
behavior  can  be  used  to  develop  an  appropriate  reduced- 
order  model  for  the  wells  —  this  model  can  be  implemented 
in  a  system  model  for  the  PCR  system.  The  implementation 
of  the  system  model  is  the  focus  of  our  current  research. 

Conclusions 

Detailed  analyses  of  PCR  well  configurations  were 
conducted  to  predict  the  contamination  between  neighboring 
wells  during  a  conventional  PCR  cycle.  Temperature 
dependent  effects  from  the  thermal  cycling  were 
incorporated  using  appropriate  diffusion  models.  The 
contamination  time  was  computed  as  the  time  required  for  5 
percent  sample  travel  between  neighboring  wells.  The 
effects  of  the  cross  sectional  geometry  as  well  as  the 
diffusion  constants  of  the  various  components  in  the  mixture 
were  analyzed  through  simulation.  The  parametric  results 
yielded  curves  that  describe  the  contamination  between 
wells  in  terms  of  the  various  parameters  that  affect  it.  These 
curves  can  then  be  applied  in  developing  optimal  designs  for 
the  wells,  as  well  as  in  the  development  of  a  model 
describing  the  behavior  of  the  PCR  well  as  a  function  of  the 
various  parameters  affecting  it.  This  model  can  then  be 
implemented  in  a  system  model  for  the  entire  device,  and  is 
the  focus  of  our  current  work. 
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Abstract 

Experimental  and  numerical  analyses  of 
multicomponent  sample  transport  are  presented  in  this 
paper,  A  novel  approach  to  confining  samples  of 
varying  electrophoretic  mobilities  during  transport  is 
demonstrated,  that  uses  a  plug  of  high  conductivity  salt 
to  reduce  the  local  field  and  thus  the  local 
electrophoretic  transport.  This  plug  confined  transport 
(PCT)  mechanism  prevents  the  separation  of  analyte 
mixtures  during  transport.  Experimental  measurements 
and  simulation  results  are  presented  that  show  good 
agreement. 

Introduction 

Electrokinetic  microfluidic  microsystems  are 
powerful  analytical  tools  for  many  applications,  such  as 
nucleic  acid  analysis,  enzyme  assays,  and 
immunoassays  [1-5].  Such  systems  have  gained 
considerable  importance  as  components  in  micron- 
scale  integrated  chemical/biochemical  analysis  or 
synthesis  systems,  also  referred  to  as  lab-on-a-chip. 
The  basic  “unit  process”  operations  in  these  systems 
are  sample  injection,  mixing,  chemical  reaction  or 
modification,  separation,  and  detection.  Assembling  a 
system  of  many  “unit  process”  nodes  requires  one  or 
more  transport  mechanisms  to  move  sample  and 
reagents  through  the  “wires”  of  the  system.  Many  of 
these  systems  rely  on  electrokinetic  physics  as  their 
transport  mechanism,  although  pressure  and  pneumatic 
applications  have  also  been  demonstrated.  Complicated 
relationships  exist  between  the  microchannel 
geometries,  the  conditions  under  which  the  devices 
operate,  and  the  behavior  of  the  multi-component  fluids 
transported  in  these  channels.  In  the  past  researchers 
have  been  forced  to  use  costly  trial  and  error  methods 
to  understand  and  design  such  microfluidic  systems. 

CAD  tools  can  be  a  valuable  aid  in  the  design  of 
microfluidic  systems.  Numerical  analyses  provide 
significant  insight  into  the  fluid  mechanics  in  these 
systems.  They  allow  the  extraction  of  material  and 
flow'  properties  that  are  generally  not  well  documented, 
or  that  vary  from  application  to  application  or  from  one 
manufacturing  technology  to  another.  Furthermore 
such  tools  help  the  designer  to  explore  a  much  larger 
space  of  designs  than  is  easily  available  from 


experiment,  and  do  so  in  a  quantitative  way  which 
enables  the  extraction  of  key  parameters  for  improved 
or  optimal  operation  of  common  microchemical  system 
components. 

Several  researchers  have  reported  CAD-based 
analyses  of  microfluidic  components  [6-9].  These 
include  components  used  in  injection  [6-7],  transport 
[8-9],  as  well  as  mixing  and  reactions.  These  analyses 
are  2D  steady-state  analyses  with  fixed  field  boundary 
conditions.  Additionally  the  numerical  analyses  were 
conducted  for  a  single  species  -  furthermore  the 
species  was  assumed  not  to  affect  the  physical  or 
electrical  properties  of  the  carrier  during  transport.  We 
have  developed  new  tools  to  enable  the  removal  of 
these  constraints  on  the  numerical  analyses.  We  have 
presented  results  for  transient  field  boundary  conditions 
showing  a  switched  injector  in  a  previous  paper  [10].  In 
this  paper,  we  extend  this  work  to  analyze 
multicomponent  mixtures,  comprising  of  species  with 
varying  electrophoretic  mobilities.  We  will  also  discuss 
the  treatment  of  high  conductivity  samples,  where  the 
sample  conductivity  significantly  affects  the  field 
within  the  carrier. 

Potential  Charged  Species 


Figure  1  :  Schematic  of  Plug  Confined  Transport  Flow, 

Although  electrokinetic  fluid  transport  has  many 
significant  advantages  such  as  allowing  simultaneous 
control  of  many  wells,  it  has  the  potential  disadvantage 
of  unwanted  separation  of  components  during 
transport.  In  this  paper  we  demonstrate  a  novel 
approach  to  confining  samples  of  varying 
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Figure  2  :  On-chip  analysis  system  showing 
network  of  interconnected  channels. 

electrophoretic  mobilities  during  transport,  by  using  a 
plug  of  high  conductivity  salt  to  reduce  the  local  field 
and  thus  the  local  electrophoretic  transport.  This  plug 
confined  transport  (PCT)  mechanism  prevents 
separation  of  the  analyte  mixtures  in  circumstances 
where  it  is  desired  to  keep  all  components  together, 
independent  of  charge.  Both  experimental  and 
modeling  data  will  be  shown  to  illustrate  this  concept. 

Objective 

There  are  two  motivations  for  this  research  in 
multicomponent  transport.  First,  is  to  aid  in  the 
development  of  microchemical  components  that 
require  mixing  or  dilution  during  electrokinetic 
transport.  Second,  is  to  aid  in  the  development  of  plug 
confined  transport  (PCT)  systems  that  use  a  high 
conductivity  salt  plug  to  transport  several  analytes  of 
different  electrophoretic  mobility.  This  novel  approach 
has  been  developed  at  Caliper  to  prevent  the 
broadening  of  analyte  mixtures  during  transport 
between  specific  sites.  A  schematic  of  the  operation  of 
PCT  is  shown  in  Figure  1.  The  plug  locally  increases 
the  conductivity  of  the  medium  and  therefore  reduces 
the  field  strength.  The  reduced  field  decreases  the 
electrophoretic  velocity  of  the  species  within  the  plug, 
which  in  turn  slows  the  migration  and  separation  of 
each  species,  confining  the  mixed  species  band  within 
the  plug  and  reducing  its  broadening. 

Experimental  Setup  and  Measurements 

Figure  2  shows  a  typical  glass  microchip  used  in  the 
validating  experiments.  Measurements  were  conducted 
in  a  straight  section  approximately  30  mm  in  length. 
The  high  conductivity  salt  plug  consists  of  a  90  mM 
HEPES  buffer  with  21mM  NaCl  at  a  pH  of  7.  The 
carrier  is  a  10  mM  HEPES  buffer  also  at  a  pH  of  7. 
Fluorescein  and  Rhodamine  were  used  as  the  markers 
to  observe  the  motion  of  the  analytes  within  the  high 


salt.  There  is  a  guard  band  of  high  conductivity  salt 
with  no  dye  around  the  high  conductivity  salt  with  dye. 
Both  the  rhodamine  and  the  fluorescein  are  carried 
along  with  the  electroosmotic  flow  of  the  carrier.  In 
addition,  the  fluorescein  also  experiences 
electrophoretic  flow  due  to  its  charge.  The 
electrophoretic  flow  is  in  the  opposite  direction  to  the 
electroosmotic  flow. 

The  fluorescent  signal  from  the  dyes  is  collected  at 
three  discrete  locations  along  the  channel  -  1.1  mm, 
13.6  mm  and  25.8  mm  from  the  plug  creation  point. 

Numerical  Formulation 

The  basic  equations  describing  the  fluid  motion  arc 
the  Navier-Stokes  equations  with  appropriate 
electromigratory  flux  terms  to  represent  the  effect  of 
the  applied  electric  field  on  the  carrier  (electroosmosis) 
and/or  the  charged  species  (electrophoresis).  The 
modeling  of  electrokinetic  effects  is  incorporated  into 
the  FlumeCAD  system.  FlumeCAD  is  an  integrated 
design  environment  consisting  of  3D  design,  modeling 
and  simulation  software  tools,  which  enable  the 
creation  and  analysis  of  complex  microfluidic  devices. 
Inherent  in  the  design  flow  implemented  in  FlumeCAD 
is  the  ability  to  translate  from  a  layout  and  process 
view  of  the  device  to  a  solid  model  and  to  continue  to  a 
3D  device  model  allowing  simulations  that  characterize 
the  various  physical  phenomena  present  in  the  device. 
The  numerical  solution  uses  a  finite  element  engine  as 
the  back-end  solver  for  the  analyses. 

Results 

In  this  section  both  experimental  and  numerical 
results  are  presented.  We  begin  by  presenting  results 
for  sample  mixing  and  dilution  to  demonstrate  the 
modeling  of  multicomponent  physics.  Results  are  then 
presented  for  separation  of  two  components  under 
electrophoresis  -  followed  finally  by  the  plug  confined 
transport  problem  -  both  experimental  and  numerical 
results  are  presented  and  compared  with  each  other. 
The  effect  of  the  guard  band  is  then  analyzed  to 
demonstrate  its  efficiency  in  confining  the  band  during 
transport. 

Mixing:  The  simplest  mixing  component  is  a  Y 
shaped  channel  as  shown  in  Figure  3.  This  component 
can  be  used  for  mixing  two  analytes  brought  in  from 
each  arm  of  the  Y  -  it  can  also  be  used  for  diluting  an 
analyte  brought  in  from  one  arm  with  pure  carrier. 
Figure  3  shows  experimental  images  for  the  latter  case 
for  electroosmotic  flow  of  a  neutral  marker  dye 
(Rhodamine).  Control  of  the  field  between  the  inlet 
ports  allows  the  varying  the  fraction  of  the  analyte  in 
the  downstream  mixture  -  and  consequently  the 
dilution  of  the  analyte.  Corresponding  simulations  for 
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Figure  3  :  Experimental  (top)  and  Simulation 
(bottom)  images  of  Mixing  of  Analyte  and  Diluent  in 
Y-channel. 


similar  cases  are  also  presented  in  Fig.  4  -  showing  the 
computed  mass  fraction  of  the  analyte.  The 
experiments  and  simulations  show  good  qualitative 
agreement  with  each  other,  as  the  figure  indicates. 

Separation:  The  separation  of  two  components 
under  the  applied  field  is  a  special  case  of  the  plug 
confined  transport  problem  -  with  no  high  conductivity 
salt  band.  Fluorescein  and  rhodamine  were  used  as  the 
marker  dyes  in  the  experiment.  Both  dyes  are 
transported  by  the  electroosmotic  flow  in  the  carrier  - 
in  addition,  fluorescein,  which  is  charged,  also 
experiences  electrophoretic  flow.  As  a  result  the  two 
species  separate.  Experimental  and  simulation  results 
are  compared  in  Figure  4  at  two  locations  along  the 
channel  -  1.1  mm  and  13.6  mm  from  the  injection 
point.  The  simulations  demonstrate  good  agreement 
with  the  experimental  measurement  -  both  in  the 
position  of  the  band  as  well  as  its  width.  The  total 
integral  under  the  band  is  not  very  well  predicted, 


especially  for  the  rhodamine  -  this  is  probably  a  result 
of  the  distortions  in  the  optical  system  from  point  to 
point.  The  good  agreement  gives  us  confidence  in  our 
ability  to  predict  the  combined  electrokinetic  transport 
of  several  analytes  -  consequently  we  now  turn  to  the 
analysis  of  plug  confined  transport  (PCT). 

Plug  confined  transport  (PCT):  Results  for  the 
PCT  problem  are  shown  in  Figures  5-7.  Experiments 
were  conducted  in  similar  conditions  to  those  discussed 
above  -  however  the  analyte  is  now  confined  in  a  high 
concentration  salt  plug  with  a  conductivity  seven  times 
higher  than  the  carrier.  The  computed  potential  field 
along  the  channel  is  shown  in  Figure  5  at  different 
instances  in  time.  The  species  mass  fraction  of  the  high 
conductivity  salt  plug  is  also  shown  in  the  figure,  along 
the  secondary  axis.  The  presence  of  the  high 
conductivity  salt  strongly  affects  the  field  as  the  figure 
shows  -  in  the  region  of  the  plug  the  field  is 
considerably  reduced.  This  region  of  low  field  travels 
downstream  along  with  the  plug.  The  electrophoresis  of 
the  charged  species,  which  is  dependent  on  the  local 
field,  is  therefore  significantly  reduced  -  thus 
preventing  electrophoretic  migration  of  the  charged 
bands  out  of  the  plug.  Diffusion,  however,  still  exists 
and  is  then  the  primary  migration  mechanism  of  the 
analyte  out  of  the  confining  band.  The  guard  bands 
around  the  sample  region  act  as  a  diffusion  barrier  to 
stop  immediate  loss  of  sample.  Once  the  analyte  leaks 
out  of  the  guard  band  into  a  region  of  high  field,  the 
overall  analyte  band  broadening  increases  dramatically 
through  electrophoresis. 

The  simulated  concentrations  of  fluorescein  and 
rhodamine  at  one  position  along  the  channel  are  shown 
in  Figure  6.  The  figure  shows  a  slice  plane  through  the 
channel  at  three  points  in  time  as  the  species  migrate 
down  the  channel.  The  concentrations  at  a  location  13.6 
mm  from  the  injection  point  are  shown  in  Figure  7  and 
compared  with  experimental  measurements  at  the  same 


Figure  4:  Experimental  and  Numerical  Measurements 
for  Separation  of  two  species  (Rhodamine  and 
Fluorescein)  .The  dashed  lines  are  experiment  and 
solid  are  simulations. 


Axial  DUtanca 


Figure  5:  Potential  Field  and  Species 
position  along  channel  for  7:1  Conductivity 
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Figure  6:  Band  Positions  of  Rhodamine  (top)  and 
Fluorescein  (bottom)  during  plug  confined  transport 
through  channel.  Images  are  at  discrete  time  intervals. 

location.  The  qualitative  agreement  between 
experiment  and  simulation  is  quite  good  -  the 
simulations  predict  the  separation  of  the  rhodamine  and 
fluorescein  peaks,  and  the  formation  of  a  second 
trailing  fluorescein  peak  due  to  the  fraction  that 
diffuses  out  of  the  edge  of  the  plug  and  is  therefore 
slowed  further  by  electrophoresis.  The  effect  of  the 
high  conductivity  salt  can  be  clearly  seen  by  comparing 
Figs.  4  and  7  for  the  cases  with  and  without  the  salt 
plug  at  the  detection  location  13.6  mm  away  from  the 
injection  point.  In  the  case  without  the  high 
conductivity  salt  plug  the  two  dyes  show  considerable 
separation,  with  the  fluorescein  arriving  almost  40 
seconds  after  the  rhodamine.  In  contrast,  for  the  plug 
confined  case,  both  peaks  arrive  almost  simultaneously, 
showing  minimal  separation. 

Conclusions 

Experimental  and  numerical  analyses  of 
multicomponent  sample  transport  are  presented  in  this 
paper.  The  mixing  of  a  sample  analyte  and  pure  carrier 
in  a  Y-shaped  mixer  is  presented  -  varying  the  field  in 
the  inlet  arms  of  the  mixer  can  control  the  dilution  of 
the  sample  with  the  carrier.  Separation  analyses  of  two 
components  -  fluorescein  and  rhodamine  are  also 
presented  under  both  electrophoresis  and 
electroosmosis.  The  negatively  charged  fluorescein  has 
a  longer  time  of  flight  than  the  neutral  rhodamine. 
Numerical  analyses  agree  quite  well  with  experiment 
both  for  position  and  width  of  the  sample  peaks. 

A  novel  approach  to  confining  the  samples  during 
transport  is  studied.  This  mechanism,  plug  confined 
transport  (PCT),  uses  a  plug  of  high  conductivity  salt  to 
transport  analytes  of  differing  electrophoretic  mobility. 
The  high  conductivity  salt  reduces  the  local  field  and 
thus  the  local  electrophoretic  transport  -  thus 
preventing  the  analytes  from  separating.  This 
mechanism  can  be  applied  in  transporting  the  analytes 
without  separation.  Experimental  and  numerical  results 
presented  here  agree  fairly  well  in  the  position  of  the 


Figure  7:  Experiment  (left)  and  Simulation  for  7:1 
Conductivity  Ratio. 

components  as  well  as  in  the  eventual  diffusion  and 
electrophoresis  of  the  chemicals  out  of  the  plug  region. 
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